Angewandte Mathematik II – Wissenschaftliches Rechnen –

Prof. Dr. rer. nat. Christian Schr¨oder

Fachbereich Ingenieurwissenschaften und Mathematik Fachhochschule Bielefeld Stand: 25. Mai 2010

1

Wie verstanden die Alten das Naturgesetz? F¨ur sie war es eine innere ” Harmonie, sozusagen statisch und unver¨anderlich; oder es war ein Idealbild, dem nachzustreben die Natur sich bem¨uhte. F¨ur uns hat ein Gesetz nicht mehr diese Bedeutung; es ist eine unver¨anderliche Beziehung zwischen der Erscheinung von heute und der von morgen; mit einem Wort: Es ist eine Differentialgleichung“ Henri Poincar´e, Der Wert der Wissenschaft

Inhaltsverzeichnis 1

Vorwort

4

2

¨ Einfuhrung

5

3

Partielle Differentialgleichungen ¨ 3.1 Ubersicht . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Grundlagen und Klassifikationen . . . . . . . . . . . . 3.3 Lineare partielle Differentialgleichungen . . . . . . . . 3.4 Lineare partielle Differentialgleichungen 2. Ordnung . 3.4.1 Klassifikation . . . . . . . . . . . . . . . . . . 3.4.2 Anfangs- und Randbedingungen . . . . . . . . 3.4.3 Korrekt gestellte Probleme . . . . . . . . . . . 3.5 L¨osung partieller Differentialgleichungen . . . . . . . 3.5.1 L¨osungsverfahren . . . . . . . . . . . . . . . . 3.5.1.1 Anwendung des Separationsansatzes

4

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

6 6 8 9 10 10 11 12 13 13 13

Numerische L¨osung partieller Differentialgleichungen 4.1 Finite-Differenzen-Methode (FDM) . . . . . . . . . . . . . . . . . . . . . 4.2 Finite-Elemente-Methode (FEM) . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Entwicklung der Finiten Elemente Methode . . . . . . . . . . . . . 4.2.2 Problemformulierung . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Variationsformulierung - Die Methode der gewichteten Residuen . . 4.2.3.1 Galerkin-Formulierung . . . . . . . . . . . . . . . . . . 4.2.4 Gebietsdiskretisierung - Finite-Elemente-Ansatzfunktionen . . . . . 4.2.4.1 Veranschaulichung via 2-Elemente-Gebietsdiskretisierung 4.2.4.2 Finite-Elemente-Ansatzfunktionen . . . . . . . . . . . . 4.2.5 Finite-Elemente-Formulierung . . . . . . . . . . . . . . . . . . . . 4.2.5.1 Assemblierung des Finite Elemente Gleichungssystems . 4.2.6 Probleml¨osung - Bestimmung der globalen Knotenwerte Q1 , Q2 , Q3 4.2.7 Finite-Elemente-L¨osung . . . . . . . . . . . . . . . . . . . . . . . 4.3 FE Ansatzfunktionen - Formfunktionen (engl. shape functions) . . . . . . . 4.3.1 Lineare Formfunktionen - Lineare Elemente“ . . . . . . . . . . . ” 4.3.2 Quadratische Elemente . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

18 18 23 24 25 26 29 29 30 33 34 36 38 38 40 41 41

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

Angewandte Mathematik II

3

Abbildungsverzeichnis 1 2 3

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Schematische Darstellung der Neumann-Randbedingungen . . . . . . . . . . . . . . . . . . . Die Abk¨uhlung eines d¨unnen Stabes der L¨ange L. Das Anfangstemperaturprofil ist durch u(x, 0) gegeben. Die Randbedingungen sind u(x, 0) = u(L, t) = 0. . . . . . . . . . . . . . . . . . . . . Dreidimensionale Darstellung der L¨osung f¨ur die Abk¨uhlung des d¨unnen Stabes der L¨ange L = 2. Das Anfangstemperaturprofil ist vorne“ durch u(x, 0) gegeben. Die Randbedingungen ” sind u(x, 0) = u(L, t) = 0. Bzgl. der Zeit t f¨allt die L¨osung nach hinten“ exponentiell ab. . . . ” Anfangs- und Randbedingung f¨ur die L¨osung der W¨armeleitungsgleichung f¨ur einen Stahlstab. Diskretisierung des L¨osungsgebiets bzgl. der x-Koordinate. . . . . . . . . . . . . . . . . . . . Diskretisierung des L¨osungsgebiets bzgl. der t-Koordinate. . . . . . . . . . . . . . . . . . . . Simultane Darstellung der Diskretisierung des L¨osungsgebiets bzgl. der x- und t-Koordinate. . Rekursives Finite-Differenzen-Verfahren zur L¨osung des eindimensionalen W¨armeleitungsproblems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Finite-Differenzen-L¨osung im Vergleich mit der analytisch exakten. . . . . . . . . . . . . . . Finite-Differenzen-L¨osung im Vergleich mit der analytisch exakten f¨ur eine zu große Schrittweite bzgl. der Zeit. Das FDM Verfahren verh¨alt sich instabil! . . . . . . . . . . . . . . . . . Gebietsdiskretisierung bei einem a) FDM und einem b) FEM Verfahren. . . . . . . . . . . . . Darstellung einer bekannten Funktion f (x) mit Hilfe eines globalen Interpolationspolynoms und mit Hilfe von Spline-Funktionen (Polynome 3. Ordnung). . . . . . . . . . . . . . . . . . Die einfachste, nicht aber unbedingt die beste, Interpolation zwischen St¨utzstellen wird durch Geradest¨ucke (Polynome 1. Grades) realisiert. . . . . . . . . . . . . . . . . . . . . . . . . . . Veranschaulichung der FEM mit Hilfe einer 2-Elemente-Gebietsdiskretisierung zur n¨aherungsweisen L¨osung des W¨armeleitungsproblems in einer Dimension. . . . . . . . . . . . . . . . . Die Ansatzfunktionen ψi (x), hergeleitet mittels geometrischer Interpretation. . . . . . . . . . . Darstellung aller Ansatzfunktionen in einem Graphen; lokale Funktionen u1 und u2 . . . . . . Einfachster Fall f¨ur die FE-L¨osung des Beispiels: Aufteilung des Gebiets in zwei Elemente Ω1 und Ω2 mit den globalen Knoten Q1 , Q2 und Q3 . . . . . . . . . . . . . . . . . . . . . . . . . Beispiel: Lineare Interpolation mittels zweier Geradenst¨ucke . . . . . . . . . . . . . . . . . . Lineare Ansatzfunktionen und die darauf basierende Darstellung der FE-N¨aherungsl¨osung uN (x) im Element Ωe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quadratische Ansatzfunktionen im Element Ωe und Koordinatentransformation der lokalen Elementr¨ander. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ab¨angigkeit der FEM-N¨aherungsl¨osung von Typ der Formfunktion. . . . . . . . . . . . . . .

12 14

16 20 20 21 21 22 22 23 24 30 31 31 32 33 36 39 41 43 43

1

Vorwort

Das vorliegende Skript zur Lehrveranstaltung Angewandte Mathematik II “ dient erg¨anzend zur Teilnahme ” an den Vorlesungen, Seminaren und Praktika, dem Selbststudium zu Hause. Auch, wenn die einzelnen Themenkomplexe ausf¨uhrlich dargestellt werden, empfehle ich zur Vertiefung des angebotenen Stoffs das Studium ¨ zus¨atzlicher Lehr- und Ubungsb¨ ucher. Als Empfehlung seien hier die folgenden B¨ucher genannt: Lehrbuch: ’Tilo Arens, Frank Hettlich, Christian Karpfinger, Ulrich Kockelkorn, Klaus Lichtenegger, Hellmuth Stachel, ”Mathematik’, Spektrum Akademischer Verlag √ Ach ja, ich verwende im Skript durchg¨angig die Bezeichung i f¨ur die imagin¨are Einheit, also i = −1! Zur Fertigstellung dieses Skripts wurde ich von Dipl.-Phys. Thomas Englisch wesentlich unterst¨utzt, wof¨ur ich ihm sehr dankbar bin. Bielefeld, Mai 2010

Christian Schr¨oder

Angewandte Mathematik II

2

5

¨ Einfuhrung e x und eine Konstante gehen die Straße herunter. Ein Ableitungsoperator ” n¨ahert sich Ihnen aus einer Seitengasse. Die Konstante bekommt große Angst und rennt davon, aber e x bleibt stehen und erkl¨art mutig, Ich bin e x ! Ich habe ” keine Angst for einem Ableitungsoperator!“. Daraufhin geht e x zum Ableix tungsoperator und sagt grinsend, Hallo! Ich bin e !“. Der Ableitungsoperator ” antwortet, Oh! Hi, ich bin ∂/∂y!“.“ ” Anonymous

Die meisten physikalischen Gesetzm¨aßigkeiten lassen sich als (in den meisten F¨allen gekoppelte partielle) Differentialgleichungen beschreiben. Hier seien als Beispiel die Maxwell-Gleichungen genannt, die die Grundlage f¨ur die gesamte Elektrotechnik bilden, von der Gleichstromtechnik bis zur H¨ochstfrequenztechnik! Die L¨osung partieller Differentialgleichungen kann in einigen F¨allen analytisch, d.h. mit Bleistift und Papier gefunden werden, in vielen F¨allen jedoch k¨onnen nur numerische L¨osungsans¨atze per Computer realisiert werden. Diese sogenannten Computersimulationen gewinnen in der Wissenschaft zunehmend dort an Bedeutung, wo herk¨ommliche Experimente oder mathematische Berechnungen nicht mehr weiterhelfen. Experimente sind in vielen F¨allen zu aufwendig oder zu riskant und mathematische Modelle erweisen sich als zu komplex, um berechenbar zu sein. Simulationen sind, in Wechselwirkung mit Experiment und Theorie, ein dritter Weg, um zu wissenschaftlichen Erkenntnissen zu gelangen. Der Aufwand f¨ur Computersimulationen ist jedoch auch nicht unerheblich! So, dass man z.B. f¨ur Klimavorhersagen, die u¨ ber viele Jahrzehnte oder Jahrhunderte reichen sollen, auf H¨ochstleistungsrechner zur¨uckgreift, die u¨ ber viele zehntausend parallel arbeitende Prozessoren verf¨ugen. Im Rahmen dieser Vorlesung werden insbesondere die mathematischen Grundlagen f¨ur die numerische L¨osung partieller Differentialgleichungen behandelt. Der Fokus liegt hierbei auf der Finiten-Elemente-Methode (FEM). Mit Hilfe der FEM ist man heutzutage in der Lage sehr komplexe technische Probleme per Computer zu simulieren. Unterst¨utzt wird dies durch eine große Auswahl an freien und kommerziellen Softwarewerkzeugen. Im Praktikum zu dieser Veranstaltung werden die Vorlesungsinhalte mit Hilfe des kommerziellen FEMProgramms COMSOL Multiphysics vertieft.

Angewandte Mathematik II

3

6

Partielle Differentialgleichungen

¨ Dieses Kapitel gibt eine Ubersicht u¨ ber den Zoo“ partieller Differentialgleichungen, die in Natur und Technik ” Verwendung finden. Wir verwenden im Text h¨aufig das u¨ bliche K¨urzel f¨ur partielle Differentialgleichungen PDE (von engl. Partial Differential Equation).

3.1

¨ Ubersicht

Um sich ein Bild von der M¨achtigkeit“ partieller Differentialgleichungen zu machen, sei im Folgenden eine ” ¨ Ubersicht u¨ ber die wichtigsten PDE im 19. und 20. Jahrhundert gegeben. ◦ 19. Jahrhundert • Laplace-Gleichung (1810) ∆u(r) = 0

(3.1)

∆u(r) = f (r)

(3.2)

• Poisson-Gleichung (1812) • Navier-Stokes-Gleichungen (1821 und 1845) ∂u(r, t) + (u(r, t) · ∇)u(r, t) − ∆u(r, t) = ∂t div u(r, t) =

−∇p

(3.3)

0

(3.4)

• Boltzmann-Gleichung (1860) ∂ f (r, v, t) F ∂ f (r, v, t) ∂ f (r, v, t) ∂ f (r, v, t) +v· + · = ∂t ∂r m ∂v ∂t

! (3.5) St

• Maxwell-Gleichungen (1873) ∂B(r, t) + rot E(r, t) = ∂t ∂D(r, t) − rot H(r, t) = ∂t div D(r, t) = div B(r, t) =

0

(3.6)

−J(r, t)

(3.7)

ρ

(3.8)

0

(3.9)

,→ Maxwell1 hat allein auf der Basis seiner“ Gleichungen elektromagnetische Wellen als L¨osung der ” PDEs vorhersagen k¨onnen! Erst viel sp¨ater im Jahre 1886 gelang die experimentelle Best¨atigung durch Heinrich Hertz2 . ◦ 20. Jahrhundert • Einsteinsche Gravitationsfeld-Gleichungen (1917) 1 8πG Rµν − gµν R + Λgµν = 4 T µν 2 c

(3.10)

,→ Bilden das Fundament der Allgemeinen Relativit¨atstheorie und pr¨agen damit unser heutiges Verst¨andnis vom Aufbau und der Dynamik des Universums. 1 James

Clerk Maxwell, schottischer Physiker (?13. Juni 1831 in Edinburgh; † 5. November 1879 in Cambridge) Rudolf Hertz, deutscher Physiker (?22. Februar 1857 in Hamburg; †1. Januar 1894 in Bonn)

2 Heinrich

Angewandte Mathematik II

7

• Schr¨odinger-Gleichung (1925) i~

! ∂ ~2 ψ(r, t) = − ∆ + V(r, t) ψ(r, t) ∂t 2m

(3.11)

,→ Sie ist das quantenmechanische Pendant zu Newtons Bewegungsgesetz und beschreibt die Physik auf atomarer Ebene. Man kann mit Fug und Recht behaupten, dass die bedeutendsten Theorien und Modelle in Natur und Technik durch PDE beschrieben werden! Sie bilden damit eine entscheidende Nahtstelle zwischen Mathematik und Ihrer Anwendungen! Hinzu kommt, dass eine bestimmte PDE h¨aufig in ganz verschiedenen Gebieten zum Einsatz kommt. Als Beispiel sei hier die Poisson-Gleichung genannt. Die L¨osung dieser PDE liefert beispielsweise • das elektrische Potential eine K¨orpers, • die vertikale Ausdehnung einer Membran, • die station¨are Temperaturverteilung in einem Volumen, oder sogar • die Anziehungskr¨afte in einem Gravitationsfeld Warum ist das so? Der Grund hierf¨ur liegt in der Herleitung“ einer PDE (irgend woher muss sie ja kommen!). ” PDEs entstammen allgemeing¨ultigen Prinzipien, d.h. vor allem physikalischen Erhaltungss¨atzen wie der Energieerhaltung, Massenerhaltung usw.! Da diese universell gelten, kann man sich vorstellen, dass auch die daraus entwickelten mathematischen Zusammenh¨ange, sprich die PDE a¨ hnlich universell eingesetzt werden k¨onnen. PDE sind wie alles auf dieser Welt auch gewissen Modeerscheinungen unterworfen. So sind derzeit besonders beliebt: • Die Kortewig-de-Vries-Gleichung (1895) ,→ zur Beschreibung von Flachwasserwellen wie z.B. Tsunamis • Black-Scholes-Gleichung (1973) ,→ Zur Bewertung von Optionen im Finanzmarkt (mit Ber¨ucksichtigun von sog. Leerverk¨aufen etc.) • Die Burger-Gleichung → zur Berechnung von Verkehrssimulationen Nach wie vor erfreuen sich PDE großer mathematischer Beliebtheit und viele Lehrst¨uhle und Fachgruppen besch¨aftigen sich mit Ihnen. Hierbei gibt es im wesentlichen zwei Forschungszweige. Der eine Zweig besch¨aftigt sich mit • qualitativen Aussagen zur L¨osbarkeit und L¨osungsmenge, ein anderer mit • quantitativen analytischen und numerischen L¨osungsmethoden Ein interessantes Problem ist das sog. Millenium-Problem“ 3 f¨ur die Navier-Stokes-Gleichungen: ” Derjenige, der die Existenz und Eindeutigkeit f¨ur ein allgemeines Anfangsproblem basierend auf den NavierStokes-Gleichungen beweisen kann, erh¨alt 1 Million US $! Kurios ist in diesem Zusammenhang die Tatsache, dass numerische N¨ahrungsl¨osungen der Navier-Stokes-Gleichungen quasi u¨ berall verwendet werden, ohne dass deren Existenz mathematisch gesichert“ist! ” 3 www.claymath.org/millenium

(Clay Mathematics Institute Cambridge, MS, USA)

Angewandte Mathematik II

3.2

8

Grundlagen und Klassifikationen

PDE: Eine PDE ist eine Bestimmungsgleichung f¨ur eine Funktion u = u(x1 , x2 , x3 , . . . , xn ) (ein skalares Feld!4 ) auf einem Gebiet Ω ⊂ ’n mit n > 1, die die • unabh¨angigen Variablen x1 , x2 , . . . , xn , • die gesuchte Funktion u und • partielle Ableitungen von u bis zu einer bestimmten Ordnung miteinander verkn¨upft! Die allgemeine Form l¨asst sich analog zur Darstellung gew¨ohnlicher Differentialgleichungen folgendermassen darstellen:     ∂u ∂u ∂u ∂2 u ∂2 u  , F  x1 , x2 , . . . , xn , u, , ,..., , . . .  = 0 (3.12) ∂x1 ∂x2 ∂xn ∂x12 ∂x1 ∂x2 Aus Gr¨unden der besseren Schreibbarkeit verwendet man h¨aufig die Abk¨urzungen: ∂u ∂2 u ; u x1 x2 ≡ ∂x1 ∂x1 ∂x2 Gleichung (3.12) schreibt sich dann kompakter als: u x1 ≡

usw.

 F x1 , x2 , . . . , xn , u, u x1 , u x2 , . . . , u xn , u x1 x1 , u x1 x2 , . . . = 0

(3.13)

Die Ordnung der PDE wird durch die h¨ochste partielle Ableitung bestimmt. Besipiele: • Die eindimensionale, zeitabh¨angige Korteweg-de-Vries-Gleichung f¨ur eine Funktion u = u(x, t): ut + 6u · u x + u xxx = 0 ⇒ Eine PDE 3. Ordnung • Die Transport-Gleichung f¨ur u = u(~u, t): ~ =0 ut + ~a · ∇u mit dem konstanten (Geschwindigkeitsvektor ~a) ⇒ Eine PDE 1. Ordnung F¨ur alle praktischen F¨alle ben¨otigt man nicht die allgemeine Form 3.12 oder 3.13, sondern z.B. • die allgemeine Form einer PDE 1.Ordnung:  F x1 , x2 , . . . , xn , u, u x1 , u x2 , . . . , u xn = 0; • die allgemeine Form einer PDE 2.Ordnung f¨ur eine zweidimensionale Funktion u = u(x, y): F(x, y, u, u x , uy , u xx , u xy , uyy ) = 0 Im letzten Fall fehlt in der Aufz¨ahlung die gemischte Ableitung u xy . Warum? Nach dem Satz von Schwarz gilt u xy = uyx ! Satz von Schwarz: Ist eine Funktion u : ’n → ’ in einer Umgebung U von ~x0 mindestens p-mal partiell differenzierbar und sind alle p-ten Ableitungen in U zumindest noch stetig so ist in ~x0 die Ableitungsreihenfolge in allen q-ten partiellen Ableitungen mit q ≤ p unerheblich. 4 Wir verwenden hier der Einfachheit halber nur die Notation f¨ ur skalare Felder; alle Beschreibungen lassen sich problemlos auf Vektorfelder u¨ bertragen.

Angewandte Mathematik II

9

Der Satz von Schwarz gilt zwar nicht allgemein, aber das soll uns weiter nicht besch¨aftigen. Dar¨uber hinaus gibt es auch Systeme von gekoppelten PDEs. Beispielhaft seien hier genannt: • Die (eindimensionale) Telegraphengleichung (zur Beschreibung der Fortpflanzung von Telegraphensignalen auf Unterseekabeln). Das Signal wird hierbei durch die zwei abh¨angigen Gr¨oßen Spannung u(x, t) und Strom i(x, t) charakterisiert, f¨ur die gilt: u x + L · it + R · i = i x + C · ut + G · u

=

0 0

Hierbei handelt es sich um ein gekoppeltes PDE System 1. Ordnung. • Die Navier-Stokes-Gleichungen als fundamentales PDE-System zur Beschreibung von Str¨omungen (→ Seminaraufgabe) • Die Maxwell-Gleichungen . . . sp¨ater mehr dazu!

3.3

Lineare partielle Differentialgleichungen

Linear partielle Differentialgleichungen spielen eine a¨ hnlich wichtige Rolle wie lineare gew¨ohnliche Differenzialgleichungen und lassen sich analog dazu definieren: lineare PDE: Eine PDE der Form:  F x1 , x2 , . . . , xn , u, u x1 , u x2 , . . . , u xn , u x1 x1 , . . . = g (x1 , x2 , . . . , xn ) heißt linear, wenn der Ausdruck F(. . . ) bzgl. der Funktion u linear ist, d.h. wenn f¨ur alle differenzierbaren Funktionen u, v : D → ’ und Faktoren λ, µ ∈ ’ gilt:  F x1 , . . . , xn , (λu + µv), (λu + µv) x1 , . . . , (λu + µv) xn , (λu + µv) x1 x1 , . . . =   = λ · F x1 , . . . , xn , u, u x1 , . . . , u xn , u x1 x1 . . . + µ · F x1 , . . . , xn , v, v x1 , . . . , v xn , v x1 x1 . . . homogene/inhomogene PDE: Wie bei gew¨ohnlichen Differentialgleichungen spricht man von einer homogenen linearen PDE, wenn g (x1 , . . . , xn ) = 0 Anderfalls ist sie inhomogen. ◦ Beispiele: Die linearen PDE 1. Ordnung skalarwertiger Funktionen u(x1 , . . . , xn ) haben alle die Gestalt: n X

a1, j (x1 , x2 , . . . , xn ) · u x j + a0 (x1 , . . . , xn ) u (x1 , . . . , xn ) = g (x1 , . . . , xn ) .

j=1

Man kann die Koeffizienten a1, j zu einem Vektor zusammenfassen und erh¨alt:   a1,1 (x1 , x2 , . . . , xn )  a (x , x , . . . , x ) n  1,2 1 2 ~a1 =  ..  .  a1,n (x1 , x2 , . . . , xn )

     .  

(3.14)

Angewandte Mathematik II

10

Unter Verwendung des Nabla-Operators l¨asst sich damit eine kompakte Schreibweise (unter Fortlassung der Argumente) finden: ~ + a0 u = g ~a1 · ∇u (3.15) In dieser Darstellung ist auch die bekannte gew¨ohnliche DGL erster Ordnung enthalten: a1 u0 + a0 u = g. Die skalaren linearen PDE 2. Ordnung haben entsprechend die Form n X

a2,i j · u xi x j +

n X

i, j=1

a1, j u x j + a0 u = g

(3.16)

··· .. . ···

(3.17)

j=1

mit der symmetrischen Matrix (warum?)   a2,11  A2 =  ...  a2,1n

3.4

a2,n1 .. . a2,nn

    

Lineare partielle Differentialgleichungen 2. Ordnung

Die meisten PDE, die in Wissenschaft und Technik relevant sind, sind 2. Ordnung! Ist das Zufall? Vielleicht nicht! In seinem nicht unumstrittenen Buch Science from Fisher Information behauptet R. B. Frieden, dass alle bekannten physikalischen Gesetze aus dem Prinzip der Minimierung der sog. Fisher-Information abgeleitet werden k¨onnen! Hierbei treten ausschließlich Operatoren 2. Ordnung als h¨ochste Terme auf. Eine interessante Behauptung! 3.4.1

Klassifikation P Der Anteil: ni, j=1 a2,i j · u xi x j aus Gleichung (3.16) wird Hauptteil einer linearen PDE 2. Ordnung genannt. Interessanterweise l¨asst sich das L¨osungsverhalten der PDEs gem¨aß der Eigenschaften der Matrix A2 aus Gleichung (3.17) vorhersagen! Es gibt drei wesentliche Typen von Gleichungen: • Elliptische PDEs: Falls A2 (~x) f¨ur alle ~x ∈ „ ausschließlich positive oder nur negative Eigenwerte hat. • Parabolische PDEs: Falls Null einfacher Eigenwert ist und alle anderen entweder positiv oder negativ sind. • Hyperbolische PDEs: Falls genau ein negativer einfacher Eigenwert vorliegt und alle weiteren positiv sind oder umgekehrt! Andere Konstellationen haben keine spezielle Bedeutung oder Namen. Die Bezeichnungen orientieren sich an den sogenannten Kegelschnitten (s. Seminaraufgabe!) Typische Vertreter sind: I Die Laplace-Gleichung: u xx + uyy + uzz = 0 ~ · ∇, ~ dem Laplace-Operator. Kurzschreibweise: ∆u = 0 mit ∆ = ∇ ⇒ elliptische PDE II Die Wellengleichung: utt − ∆u = 0 ⇒ hyperbolische PDE III Diffusionsgleichung oder W¨armeleitungsgleichung: ut − ∆u = 0 ⇒ parabolische PDE

Angewandte Mathematik II

11

Diese Gleichungen sind Prototypen ganzer Klassen von PDEs! Jede Klasse hat ein charakteristisches L¨osungsverhalten. I Elliptische PDE: Station¨ares Verhalten { Information (d.h. die Randbedingungen!) wird ohne Zeitverlust, d.h. instantan, zu allen Punkten im Gebiet Ω transportiert. II Hyperbolische PDE: Periodisches Verhalten { Information propagiert“ wellenartig zu Punkten (nicht notwendigerweise alle!) im Gebiet Ω. ” III Parabolische PDE: Exponentielles Wachstum oder Zerfall { Information diffundiert“ zu allen Punkten im Gebiet Ω. ” 3.4.2

Anfangs- und Randbedingungen

Genau wie bei gew¨ohnlichen Differentialgleichungen ben¨otigt man auch bei PDE Anfangs- und Randbedingungen, um aus der allgemeinen L¨osung eine spezielle L¨osung zu bestimmen. Insbesondere f¨ur die Berechnung numerischer N¨aherungsl¨osungen sind Randbedingungen wichtig. Aber dazu sp¨ater mehr. Anfangsbedingungen : Bei einer Funktion u(~x, t) spricht man von Anfangsbedingungen, wenn gilt: u(~x, t0 ) = f (~x) oder

∂u (~x, t0 ) = ~g(~x) ∂t

(3.18)

an einer Stelle t0 ∈ ’, wenn t im Sinne einer Zeitabh¨angigkeit interpretiert werden kann. Beachten Sie! Die Anfangsbedingungen (3.18) sind grunds¨atzlich immer noch ortsabh¨angig! Falls f (~x) = 0 bzw. g(~x) = 0 nennt man die Anfangsbedingungen homogen. Dirichlet-Randbedingungen 5 : Man kann f¨ur u(~x, t) die Funktionswerte auf dem Rand ∂Ω vorgegeben werden: u(~x, t) = g(~x) ; ~x ∈ ∂Ω (3.19) Man bezeichnet die Randbedingungen (3.19) als Dirichlet-Randbedingungen. Beachten Sie! Die Randbedingungen (3.19) sind grunds¨atzlich immer noch zeitabh¨angig! ~ x, t), in Richtung des NorNeumann-Randbedingungen 6 : Man kann auch die Ableitung von u(~x, t), also ∇u(~ malenvektors ~n zum Rand von Ω angeben: ~ x, t) = ~u · ∇u(~

∂u(~x, t) = g(~x); ∂~n

~x ∈ ∂Ω

Diese Randbedingung nennt man Neumann-Randbedingung (siehe Abb.1). Robin-Randbedingung Dar¨uber hinaus gibt es Mischformen der Art: αu(~x) + β

∂u(x) = g(~x) ∂~n

; ~x ∈ ∂Ω

Man nennt diese Robin-Randbedingungen (wobei Robin hierbei franz¨osisch ausgesprochen wird!). 5 benannt nach dem deutschen Mathematiker Johann Peter Gustav Lejeune Dirichlet ?13. Februar 1805 in D¨ uren; †5. Mai 1859 in G¨ottingen) 6 benannt nach dem deutschen Mathematiker Carl Gottfried Neumann (?7. Mai 1832 in K¨ onigsberg (Preußen); †27. M¨arz 1925 in Leipzig)

Angewandte Mathematik II

12

u n

n

u(x)



Abbildung 1: Schematische Darstellung der Neumann-Randbedingungen

Wie im Fall der gew¨ohnlichen Differentialgleichungen haben Anfangs- und Randbedingungen entscheidenden Einfluss auf die L¨osungen bzw. die Anwendung von L¨osungsmethoden! Gibt man beispielsweise nicht” kompatible“ Randbedingungen an, so findet man unter Umst¨anden keine spezielle L¨osung, obwohl eine allgemeine L¨osung existiert! Passiert dies in Comsol, so tauchen kryptische“ Fehlermeldungen auf, die nicht ” unbedingt auf dieses Problem schließen lassen, aber dazu im Praktikum mer! Andernfalls kann eine falsch“ ” definierte Randbedingung, also eine, die zwar mathematisch korrekt, aber physikalisch unsinnig ist, zu einer speziellen L¨osung f¨uhren, die jedoch nichts mit der Wirklichkeit zu tun hat! Insgesamt nennt man eine PDE mit ihren Anfangs-, Rand-, oder Anfangs- und Randbedingungen: • Anfangswertprobleme (AWP) • Randwertprobleme (RWP), oder • Anfangs-Randwert-Probleme (ARWP) 3.4.3

Korrekt gestellte Probleme

Die Empfindlichkeit“ einer L¨osung hinsichtlich der Wahl ihrer Rand- und Anfangsbedingungen hat zum Be” griff des korrekt-gestellten Problems gef¨uhrt. Man kann sagen, dass f¨ur ein korrekt-gestelltes Problem die Anfangs- und Randbedingungen so formuliert sein m¨ussen, dass das Gesamtproblem garantiert eine eindeutige L¨osung von praktischer Relevanz besitzt! Insbesondere darf das Problem nicht sensitiv gegen¨uber kleinen ¨ Anderungen der Anfangs- und Randbedingungen sein! Definition (nach J.S. Hadamard (1865 - 1963)): Korrekt-gestelltes Problem : Ein AWP, RWP oder ARWP heißt korrekt-gestellt (engl. well-posed), wenn 1 eine lokale L¨osung des Problems existiert, 2 die lokale L¨osung eindeutig ist, ¨ 3 die lokale L¨osung stetig von den Anfangs- und Randbedingungen abh¨angt. D.h. kleine Anderungen ¨ in den Anfangs- und Randbedingungen nur kleine Anderungen in der L¨osung verursachen! Beachten Sie, dass diese Bedingungen nur lokaler Natur sind; d.h. in einem gegebenem Gebiet existiert zwar eine eindeutige lokale L¨osung, das globale Problem kann jedoch auch weitere L¨osungen besitzen! Zur Orientierung: Um ein korrekt-gestelltes Problem zu formulieren, ben¨otigt man auf jeden Fall je Zeitableitung eine Anfangsbedingung und je Ortsableitung eine Randbedingung. Beispiel: F¨ur die eindimensionale Wellengleichung utt − u xx = 0 ben¨otigt man 2 Anfangsbedingungen und 2 Randbedingungen, um ein korrekt-gestelltes Problem zu formulieren. (s. Seminaraufgabe!)

Angewandte Mathematik II

3.5

13

L¨osung partieller Differentialgleichungen

Analog zum Fall der gew¨ohnlichen Differentialgleichungen bezeichnet man als L¨osung eine Funktion, die • gen¨ugend oft partiell differenzierbar ist, • eingesetzt in die PDE ergibt sich 0 = 0 im Gebiet Ω; • man unterscheidet zwischen allgemeiner und spezieller L¨osung, wobei die allgemeine f¨ur die Praxis keine Relevanz besitzt. 3.5.1

L¨osungsverfahren

Es gibt eine große Vielfalt an analytischen L¨osungsverfahren, die f¨ur sich jedoch nicht sehr allgemein sind und nur auf bestimmte Typen von PDE anwendbar sind. Die wichtigsten sind hierbei: Substitutionsverfahren : Hierbei werden die Variablen durch eine Transformation geschickt ersetzt, so dass sich eine vereinfachte partielle oder sogar gew¨ohnliche Differentialgleichungen ergeben. Charakteristikenverfahren : Man bestimmt hierbei eine Parametrisierung des Graphen der L¨osung. Dieses Verfahren ist auf sog. quasi-lineare Differentialgleichungen erster Ordnung beschr¨ankt (vgl. Richtungsfeld einer expliziten gew¨ohnlichen DGL 1. Ordnung). Dies f¨uhrt auf ein System gew¨ohnlicher Differentialgleichungen. Potenzialtheorie : Man verwendet hier f¨ur den Ansatz einen Potenzialausdruck oder Summe von Potenzialen bei einem Randwertproblem. Dies f¨ur auf eine Integralgleichung. Variationsformulierung : Entweder direkt durch partielle Integration oder durch Formulierung als ein Optimierungsproblem kann eine Variationsformulierung eines Anfangs- oder Randwertproblems gefunden werden. Separationsansatz : Ansatz als Produkt von Funktionen, die jeweils nur von einer unabh¨angigen Variable abh¨angen. Dieser Ansatz f¨uhrt auf ein System von (ungekoppelten) gew¨ohnlichen Differentialgleichungen. Integraltransformationen : Eine Transformation bez¨uglich einer oder mehrerer Variablen mit Hilfe der Laplace- oder Fouriertransformation f¨uhrt auf eine auf eine gew¨ohnliche Differentialgleichung im Bildraum, die dort mit den bekannten Verfahren gel¨ost werden kann. Anschließend wird die L¨osung zur¨uck transformiert, womit man die L¨osung der PDE erh¨alt. Wie bei jeder anderen Anwendung von Integraltransformationen, funktioniert dieses Verfahren nur, wenn man die R¨ucktransformation auch wirklich durchf¨uhren kann. Ein Beispiel hierzu ist als Seminaraufgabe gegeben. Den Ansatz mit Hilfe von Integraltransformationen haben wir schon in Angewandte Mathematik I“ studiert. ” Der Separationsansatz ist ebenfalls bekannt (Erinnern Sie sich an das Problem der Eisversorgung Alexander ” des Großen“!). Wir wollen aber dennoch anhand eines Beispiels n¨aher darauf eingehen. 3.5.1.1 Anwendung des Separationsansatzes Gesucht sei die L¨osung zum ARWP f¨ur die Diffusionsgleichung in einer Dimension, d.h. f¨ur u = u(x, t) mit: ut − a2 u xx = 0 ;

a2 > 0 (Konstante) ,

mit: 0 ≤ x ≤ L ; t > 0

zu den Randwerten: u(0, t) = u(L, t) = 0

(→ Dirichlet-Randbedingungen)

Angewandte Mathematik II

14

Die Anfangsbedingung sei

x  π . L Dieses ARWP beschreibt beispielsweise die Abk¨uhlung eines d¨unnen Stabes der L¨ange L, der an den Enden auf konstanter Temperatur gehalten wird. Beachten Sie, dass die Rand- und Anfangsbedingungen immer konsistent sein m¨ussen! D.h. in diesem Fall, dass die Anfangsbedingung an den R¨andern auch den Wert Null haben muss. In einem Werkzeug wie Comsol kann man u¨ blicherweise problemlos inkonsistente Rand- und Anfangsbedingungen getrennt von einander eintragen - ohne, dass sich die Software beschwert! Dies f¨uhrt beim Starten der Simulation unter Umst¨anden zu kryptischen“ Fehlermeldungen, aus denen nicht ersichtlich wird, was hier schief gelaufen ist. ” u(x, 0) = sin

u(x,0) u(0,t)

u(L,t) L

0

Abbildung 2: Die Abk¨uhlung eines d¨unnen Stabes der L¨ange L. Das Anfangstemperaturprofil ist durch u(x, 0) gegeben. Die Randbedingungen sind u(x, 0) = u(L, t) = 0. F¨ur den L¨osungsansatz w¨ahlt man nun u(x, t) = v(t) · w(x). Einsetzen in die PDE mit

ut = v (t) · w(x) 0



und u xx = v(t) · w00 (x) ergibt:

v0 (t) · w(x) − a2 v(t) · w00 (x) = 0 | : (w(x)v(x)) 00 v0 (t) 2 w (x) =a (f¨ur w(x) , 0 ; v(t) , 0 !) v(t) w(x)

Idee: Die linke Seite h¨angt nur von t ab und die rechte Seite h¨angt nur von x ab! ⇒ Damit diese Bedingung erf¨ullt ist, m¨ussen beide Ausdr¨ucke gleich derselben Konstanten, z.B. a2 k sein7 Daraus folgt:

v0 (t) = a2 k v(t)

und

a2

w00 (x) w00 (x) = a2 k ⇔ = k; w(x) w(x)

k∈ƒ

⇒Man erh¨alt somit zwei separierte gew¨ohnliche Differenzialgleichungen! 2

I v0 (t) = a2 k · v(t) ⇒ Allgemeine L¨osung: v(t) = C1 · ea kt (gew¨ohnliche Differentialgleichunge 1. Ordnung, l¨osbar mit dem Verfahren Trennung der Variablen“) ” √



II w00 (x) − kw(x) = 0 ⇒ Allgemeine L¨osung: w(x) = C2 · e kx + C3 · e− kx (gew¨ohnliche lineare, homogene Differenzialgleichung 2. Ordnung mit konstanten Koeffizienten) (⇒ Beachten Sie, man erh¨alt drei Integrationskonstanten, d.h. f¨ur jede partielle Ableitung eine!) Einsetzen der Randbedingunen: I u(0, t) = v(t) · w(0) = 0 II u(L, t) = v(t) · w(L) = 0 7 Man h¨ atte auch einfach nur k nehmen k¨onnen; die Verwendung von a2 k f¨uhrt schlicht nur zu einer einfacheren Schreibweise! Probieren Sie es aus!

Angewandte Mathematik II

15

Annahme: v(t) , 0 f¨ur t > 0 (f¨ur t = 0 gibt es ja die Anfangsbedingung!) { andernfalls erg¨abe sich sofort die L¨osung: u(x, t) = 0 f¨ur t > 0 ⇒ Diese L¨osung bezeichnet man als triviale L¨osung“, da sie zwar existiert, aber im Sinne der Problemstellung v¨ollig uninteressant ist! ” Damit ergibt sich zwangsl¨aufig w(0) = 0 und w(L) = 0. Setzt man dies nun in die allgemeinen L¨osungen ein, so ergibt sich I:

= C2 + C3 = 0

w(0)

II: w(L) = C2 · e ⇒ C3

=



0

=



0

=



kL

und

+ C3 · e−

√ kL

=0

−C2  √ √  C2 e kL − e− kL √  2C2 sinh kL

√ • F¨ur k > 0 muss C2 = 0 , da der sinh( kL) , 0. Daraus folgt C3 = 0 • F¨ur k < 0 ist w(x) = C2 ei



w(x) ≡ 0



∀x



und wir erhalten wieder eine triviale L¨osung!



+ C3 e−i kx (wobei k jetzt nur noch den Betrag von k, als |k| sei) √  √  ¨ e3 cos kx e2 sin kx + C (Ubergang zur reellen Basis) w(x) = C



kx

Einsetzen der Randbedingungen: =

0



II: w(L) =

0



I:

w(0)

e3 = 0, da cos (0) = 1 C √  e2 sin k · L 0=C

e2 , 0 ist Gleichung (II) immer dann erf¨ullt, wenn: F¨ur ein C √ kL = n · π f¨ur n ∈ Ž bzw. ⇒



k=

und damit

n·π L

(?)

F¨ur jedes n ∈ Ž gibt es eine L¨osung: wn (x) = Cn · sin

 nπ  x L

mit Cn ∈ ’

Die hierzu passenden“ zeitlichen L¨osungen ergeben sich mit v(t) = C1 e−a ” Bedingung (?):  nπ 2 √ nπ k = ⇔ k= L L 2 ·t −( nπa ) mit Bn ∈ ’ ⇒ vn (t) = Bn · e L

2

kt

durch Ersetzen von k mit der

S¨amtliche L¨osungen der PDE ergeben sich nun zu un (x, t) = =

vn (t) · wn (x)  nπ  nπa 2 Dn · e−( L ) t sin x ; mit Dn = Bn · Cn , n ∈ Ž. L

Da die PDE linear ist, ist auch die Linearkombinationen aller un (x, t) eine L¨osung der PDE: ⇒

u(x, t) =

∞ X n=1

Dn · e−(

nπa 2 t L

)

sin

 nπ  x , L

Angewandte Mathematik II

16

d.h. die L¨osung u(x, t) ist schlicht eine Fourier-Reihe bzgl. x! Setzt man die Anfangsbedingung: u(x, 0) = sin Lx π ein, so erh¨alt man  nπ  x  ! un (x, 0) = Dn sin x = sin π , (3.20) L L woraus folgt, dass n = 1 und damit D1 = 1 ist. Man erh¨alt somit π  πa 2 u(x, t) = e−( L ) t sin x L als L¨osung des gegebenen Anfangsrandwertproblems.

Abbildung 3: Dreidimensionale Darstellung der L¨osung f¨ur die Abk¨uhlung des d¨unnen Stabes der L¨ange L = 2. Das Anfangstemperaturprofil ist vorne“ durch u(x, 0) gegeben. Die Randbedingungen sind u(x, 0) = u(L, t) = ” 0. Bzgl. der Zeit t f¨allt die L¨osung nach hinten“ exponentiell ab. ” Wir wollen noch einmal auf die Ausgangssituation zur¨uck kommen. Dieses ARWP beschreibt die Abk¨uhlung eines d¨unnen Stabes der L¨ange L, der an den Enden auf konstanter Temperatur gehalten wird. Was heißt in diesem Zusammenhang d¨unn“? Genau solche Fragen sind zentraler Bestandteil bei der mathematischen Mo” dellierung realer Situationen. Und es sind insbesondere die Randbedingungen, die die Grenzen eines Modells hinsichtlich der Anwendbarkeit festlegen (denken Sie an die Aufgabe mit den drei Lampen und Schaltern!)! In diesem eindimensionalen Fall gibt es lediglich die beiden Enden des Stabes, die als R¨ander Werte tragen k¨onnen. Nur hier kann ein Informationsaustausch (z.B. konstante Temperatur oder ein W¨armefluss) mit der Umgebung stattfinden. Die spannende Frage ist nun wie sich dies mit einer realen Situation vereinbaren l¨asst? Stellt man sich einen sehr d¨unnen Metalldraht vor, der auf einer Seite erhitzt wird, so ist dies mitnichten die Situation, die hier mathematisch modelliert wird! Warum nicht? Der Metalldraht wird, und sei er noch so d¨unn, nach einiger Zeit insgesamt warm werden, d.h. man kann ihn anfassen und m¨oglicherweise schmerzhaft f¨uhlen wie die W¨arme durch den Draht transportiert wird. Dies ist ein Informationsaustausch u¨ ber die Staboberfl¨ache und widerspricht vollkommen der Modellannahme, nach der wir gar keine Staboberfl¨ache haben d¨urften! Wir haben hier folgendes Problem. Der Begriff d¨unn bezieht sich nur auf die Formulierung der Diffusionsgleichung in einer Dimension und darf nicht unbedingt direkt in die Realit¨at u¨ bertragen werden! Denn in der Realit¨at hat selbst ein d¨unner Draht eine Oberfl¨ache und ist damit dreidimensional! Wenn wir nun unser sch¨ones, analytisch l¨osbares Modell retten wollen, dann m¨ussen wir uns eine passende reale Situation konstruieren. Die Nichtexistenz einer Oberfl¨ache im Modell bedeutet, dass hier kein Informationsaustausch mit der Umgebung stattfinden kann, d.h. der Stab hinsichtlich seiner Oberfl¨ache (aber nicht an den Enden!) isoliert ist. Bezogen auf die Realit¨at ist dieses Modell also nur dann anwendbar, wenn der Stab entsprechend isoliert ist, z.B. durch einen

Angewandte Mathematik II

17

Styropor-Block verl¨auft, wobei streng genommen auch dies nur eine N¨aherung der mathematisch geforderten Isolation darstellt. Wenn nun anders herum die physikalische Situation eines Stabes, der in Styropor eingebettet ist, vorliegt, bedeutet dies nun zwangsl¨aufig, dass wir diese Situation mit Hilfe der eindimensionalen Diffusionsgleichung beschreiben m¨ussen? Nein, ganz und gar nicht. Wir k¨onnen direkt ein dreidimensionales Modell verwenden (das sich dann aber nur noch per Computer l¨osen l¨asst). Allerdings m¨ussten wir dann hier explizit die Isolierung der Oberfl¨ache (bis auf die Enden) als Randbedingung angeben. Was w¨are damit gewonnen? Nichts, im Gegenteil! Ein dreidimenionales Modell schl¨agt sich in einem gewaltig h¨oherem Rechenaufwand nieder - und das mit dem selben Ergebnis! Sobald n¨amlich die Oberfl¨ache als isolierend betrachtet wird, spielen W¨armefl¨usse senkrecht zur Oberfl¨ache keine Rolle mehr und damit reduziert sich der W¨armetransport auf eine Dimension, so dass ein eindimensionales Modell v¨ollig ausreicht! Die Modellierung der Randbedingungen ist also ganz entscheidend f¨ur die Verkn¨upfung von Modell und Realit¨at und erfordert besondere Aufmerksamkeit! Wir wollen dies im folgenden Kapitel n¨aher beleuchten.

Angewandte Mathematik II

4

18

Numerische L¨osung partieller Differentialgleichungen

Einige Warnungen vorweg! Im Gegensatz zu analytischen L¨osungen liefern numerische L¨osungsverfahren keine Formel bzw. Gleichung f¨ur die L¨osung, sondern nur N¨aherungen f¨ur die Funktionswerte im betrachteten Gebiet. Dies hat ganz entscheidende Konsequenzen. So kann man beispielsweise mit einer analytischen L¨osung (siehe Beispiel aus dem vorherigen Kapitel) den Zustand eines Systems zu jeder Zeit (sofern zeitabh¨angig) und an jedem Ort direkt durch Einsetzen von Zahlen f¨ur die unabh¨angigen Variablen ermitteln. Bei einer numerischen L¨osung geht das nicht! Der Aufwand ist hier viel h¨oher! Man muss sich hier der gew¨unschten L¨osung im wahrsten Sinn des Wortes n¨ahern, d.h. man muss z.B. f¨ur einen bestimmten Zeitpunkt eine Menge von vorherigen Zeitpunkten bestimmen oder f¨ur einen bestimmten Ort auch alle Nachbarorte bzw. das vollst¨andige Gebiet berechnet haben. Hinzu kommt, dass numerische Verfahren (fast) immer eine L¨osung“liefern, die ” jedoch beliebig falsch sein kann! Gr¨unde daf¨ur sind: • Rundungsfehler (man spricht hier vom Maschinenproblem“) ” • Fehlerhafte N¨aherung exakter mathematischer Operatoren (sog. Methodenfehler“) ” ¨ • Abbildungsfehler“ bei der Uberf¨ uhrung der Realit¨at in ein Modell. Allein bzgl. der Geometrie passieren ” hier Fehler aufgrund der Diskretisierung des betrachteten Gebiets besonders an den R¨andern.) Der Aufwand f¨ur die Bestimmung einer numerischen L¨osung steigt mit der Komplexit¨at des Problems/Modells (nur selten linear!) an, so dass es passieren kann, dass Rechnungen praktisch undurchf¨uhrbar werden! Hier st¨oßt man auf ein generelles Dilemma. Je wirklichkeitsn¨aher ein Modell werden soll, desto aufwendiger wird seine numerische L¨osung. Man muss also in jedem Fall Abstriche machen, sonst kann es schnell passieren, dass eine Modellrechnung v¨ollig unn¨utz ist (z.B. wenn die Berechnung der Wettervorhersage f¨ur den n¨achsten Tag eine Woche dauert ...) Es gibt eine Vielzahl numerischer Verfahren, die sich in zwei Kategorien einteilen lassen: • Direkte Verfahren, d.h. die direkte Verwendung der PDE zur L¨osung des Problems; hierbei werden kontinuierlicher Ausdr¨ucke, insbesondere Differentialoperatoren durch diskrete N¨aherungen auf einem Gitter ersetzt. Wichtigster Vertreter ist hier das Finite-Differenzen-Verfahren (FDM) • Projektionsverfahren. Es findet eine Umformulierung bzw. Projektion des Problems bzw. der PDE statt, um zu einer einfacheren Beschreibung zu kommen. Hier ist die Finite-Elemente-Methode der wichtigste Vertreter. Wir wollen die beiden genannten Verfahren im Folgenden vorstellen, wobei sich die Diskussion der FDM in Grenzen h¨alt. Unser Fokus wird auf der FEM liegen, die die mathematische Grundlage f¨ur das eingesetzte Werkzeug COMSOL Multiphysics darstellt.

4.1

Finite-Differenzen-Methode (FDM)

F¨ur die FDM wird zun¨achst die Definitionsmenge des Problems wird mit einem (1D, 2D, 3D) kartesischen Gitter u¨ berzogen, d.h. diskretisiert. Auf diesem Gitter werden Ersatzausdr¨ucke“ f¨ur die exakten Differential” operatoren definiert. Die differentiellen Ausdr¨ucke werden durch endliche“ Differenzenquotienten (daher der ” Ausdruck finite Differenzen“) bzgl. der unabh¨angigen Variablen ersetzt. Prinzipiell kann man hierbei den Feh” ler reduzieren, in dem man immer kompliziertere Differenzenquotienten definiert. Allerdings stehen auch hier Aufwand und Machbarkeit in Konkurrenz. Aus der PDE wird ein (endlich dimensionales) Gleichungssystem, das die Funktionswerte der L¨osung an den Gitterpunkten als unbekannte Variablen enth¨alt. Man kann schon ahnen, dass die Genauigkeit der Methode mit der Dichtheit“ des Gitters zusammen h¨angt. Je dichter das Gitter, ” desto gr¨oßer das zu l¨osende Gleichungssystem. Da die Dimension des Problems als Potenz in die Anzahl der Unbekannten eingeht (z.B.: 1D:10 Gitterpunkte ⇒ 3D: 103 = 1000 Gitterpunkte!) kommt hier sehr schnell

Angewandte Mathematik II

19

zu sehr großen Gleichungssystemen. Konstruktionsbedingt (die N¨aherungen der Differentialoperatoren sind lokal, d.h. es werden nur benachbarte Punkte einbezogen) sind die FDM-Gleichungssysteme jedoch d¨unn besetzt, d.h. viele Eintr¨age sind Null. F¨ur solche Gleichungssysteme gibt es eine Vielzahl sehr effizienter exakter und (f¨ur sehr große Systeme) n¨aherungsweiser L¨osungsverfahren. Wir wollen im Folgenden das Verfahren am beliebten Beispiel der L¨osung der W¨armeleitungsgleichung in einer Dimension erl¨autern. Sei u(x, t) die Temperatur eines K¨orpers an einem Ort x zur Zeit t, so beschreibt k · ∆u(x, t) = u(x, t)t jegliche Ver¨anderung von u(x, t) in Raum und Zeit. Zur Erinnerung, die 2. Ortsableitung (d.h. der LaplaceOperator) ist definiert als ∂2 u(x, t) ∂2 u(x, t) ∂2 u(x, t) + + , ∆u(x, t) = ∂x2 ∂y2 ∂z2 die vorkommende 1. Zeitableitung lautet ausf¨uhrlich geschrieben u(x, t)t =

∂u(x, t) . ∂t

Der Wert k ist die konstante Temperaturleitf¨ahigkeit und damit die einzige(!) Materialgr¨oße in diesem Problem (die Temperaturleitf¨ahigkeit f¨ur Stahl ist k = 0.15cm2 /s). Wie im letzten Kapitel wollen wir der einfachheithalber die L¨osung der W¨armeleitungsgleichung in einer Dimension betrachten. Also auch hier gehen wir von einem d¨unnen Stahlstab der L¨ange L aus, der bis auf die Enden isoliert sei (also nur ein W¨armefluss in x-Richtung m¨oglich ist!). Das Problem reduziert sich damit auf die L¨osung von k

∂2 u ∂u = , ∂t ∂x2

mit der Anfangsbedingung u(x, t) t=0 = f (x)

f¨ur alle Orte x ∈ Ω

und der Dirichlet-Randbedingungen u(x, t) x=0 u(x, t) x=L

=

0

=

0

und

Wir wollen nun die Approximation der W¨armeleitungsgleichung mit Hilfe der FDM betrachten. Alle Varianten des Finite-Diffenzen-Verfahrens approximieren eine (exakte) L¨osung u(x, t) auf einem vorgegebenen Gebiet Ω an diskreten St¨utzstellen! Die Vorgehensweise ist folgendermaßen: ¨ • Uber das Gebiet Ω wird ein Gitternetz“ gespannt. ” ¨ • Alle in der PDE auftretenden Funktionswerte und Ableitungen werden durch finite Ausdrucke ersetzt. Wie in Abb. 5 gezeigt f¨uhrt dies zu einer Approximation der Funktionswerte an diskreten St¨utzstellen gem¨aß u(x, t) ≈ u(xi , t) = ui . Die 1. und 2. Ortsableitung l¨aßt sich mit Hilfe der folgenden Differenzenquotienten darstellen: ∂u(xi , t) ∂x ∂u(xi , t) ∂x 2 ∂ u(xi , t) ∂x2

≈ ≈ ≈

1 [ui+1 − ui ] (Vorw¨artsdifferenz) hx 1 [ui − ui−1 ] (R¨uckw¨artssdifferenz) hx 1 [ui+1 − 2ui + ui−1 ] (Zentrale Differenz) h2x

Angewandte Mathematik II

20

Abbildung 4: Anfangs- und Randbedingung f¨ur die L¨osung der W¨armeleitungsgleichung f¨ur einen Stahlstab.

Abbildung 5: Diskretisierung des L¨osungsgebiets bzgl. der x-Koordinate.

F¨ur die Approximation der Zeitabh¨angigkeit gilt analog (siehe Abb. 6) f¨ur die Funktionswerte u(x, t) ≈ u(x, t j ) = u j sowie f¨ur die 1. Zeitableitung ∂u(x, t j ) ∂t ∂u(x, t j ) ∂t

≈ ≈

1 j+1 [u − u j ] (Vorw¨artsdifferenz), ht 1 j [u − u j−1 ] (R¨uckw¨artssdifferenz). ht

Man kann die Approximation der Orts- und Zeitabh¨angigkeit simulatan in einem im Raum-Zeit-Gitter“ dar” stellen (siehe Abb. 7). Mit 1 ∂2 u(xi , t) ≈ 2 [ui+1 − 2ui + ui−1 ] und 2 ∂x hx ∂u(x, t j ) 1 ≈ [u j+1 − u j ] ∂t ht kann man die W¨armeleitungsgleichung k

∂2 u ∂u = ∂t ∂x2

in die Finite-Differenzen-Form umschreiben: k

1 j 1 j (u − 2uij + ui−1 ) = (uij+1 − uij ) ht h2x i+1

Angewandte Mathematik II

21

Abbildung 6: Diskretisierung des L¨osungsgebiets bzgl. der t-Koordinate.

Abbildung 7: Simultane Darstellung der Diskretisierung des L¨osungsgebiets bzgl. der x- und t-Koordinate.

Aufl¨osen dieser Gleichhung nach uij+1 ergibt ein sog. explizites Differenzenverfahren j j uij+1 = r(ui+1 + ui−1 ) + (1 − 2r)uij

mit

r=k

ht , h2x

das die Berechnung der unbekannten N¨aherungswerte mit Hilfe bereits bekannter (vorher berechneter) an den St¨utzstellen erlaubt. In Abb. 8 ist das Berechnungsschema gemeinsam mit den Anfangs- und Randbedingungen (als blaue Punkte) • Anfangsbedingung: u1i = f (xi ), i = 1, 2, .., N • Randbedingung links: u1j = 0, j = 1, 2, .., M • Randbedingung rechts: uNj = 0, j = 1, 2, .., M dargestellt. Man erkennt sofort, dass sich hinter Gleichung (4.1) ein rekursives Berechnungsverfahren verbirgt. Die L¨osung f¨ur das gesamte Gebiet wird zeilenweise von links nach rechts und von unten nach oben berechnet. Verwendet man beispielsweise h x = 0.25 und r = 0.5, so ergibt sich f¨ur die Schrittweite bzgl. der Zeit ht = 0.206. Damit kann man die approximative L¨osung bestimmen (siehe Abb. 9). Man erkennt im Vergleich mit der analytisch exakten L¨osung, dass die FDM eine brauchbare N¨aherung liefert. Man ist nun versucht, z.B. die eine gr¨oßere Schrittweite bzgl. der Zeit zu w¨ahlen, um Zeit bei der numerischen Berechnung zu sparen. W¨ahlt man beispielsweise r = 0.6, so ergibt sich ht = 0.248. F¨ur diese Parameter liefert das vorgestellte Verfahren keine verl¨aßlichen Werte mehr wie in Abb. 10 zu sehen ist! Man kann zeigen, dass das Verfahren nur unter der Bedingung r ≤ 0.5 bzw. ht ≤ k/2 h2x stabile Ergebnisse produziert, dem entsprechned ist es f¨ur r > 0.5 bzw. ht > k/2 h2x instabil!

Angewandte Mathematik II

22

Abbildung 8: Rekursives Finite-Differenzen-Verfahren zur L¨osung des eindimensionalen W¨armeleitungsproblems.

Abbildung 9: Finite-Differenzen-L¨osung im Vergleich mit der analytisch exakten.

Das Problem der Instabilit¨at l¨aßt sich mit einem einfachen Trick beheben! Anstatt in der Finite-DifferenzenForm 1 j 1 j+1 j k 2 (ui+1 − 2uij + ui−1 ) = (ui−1 − uij ) ht hx die linke Seite nur zur Zeit t (d.h. Gitterlinie j) zu betrachten, bildet man den Mittelwert aus den Ausdr¨ucken zur Zeit t und t + 1 (d.h. j und j + 1) und erh¨alt k

i 1 j+1 1 h j j j+1 j+1 − uij ) (ui+1 − 2uij + ui−1 ) + (ui+1 − 2uij+1 + ui−1 ) = (ui−1 2 ht 2h x

Aufl¨osen der Gleichung nach den unbekannten u j+1 f¨uhrt auf ein implizites Differenzenverfahren (das sog. Crank-Nicolson-Verfahren). j+1 j+1 j j −rui+1 + 2(1 + r)uij+1 − rui−1 = rui+1 + 2(1 − r)uij + rui+1

Man kann diese Gleichung mit Hilfe von Matrizen kompakter darstellen als Cu j+1 = Au j

Angewandte Mathematik II

23

Abbildung 10: Finite-Differenzen-L¨osung im Vergleich mit der analytisch exakten f¨ur eine zu große Schrittweite bzgl. der Zeit. Das FDM Verfahren verh¨alt sich instabil!

mit

u j+1

 j+1    u1   1  j+1  −r  u2    0  u j+1    3  =  .  ; C =  .  ..  ..      0 u j+1     N−1 j+1  0 uN

0 0 0 2(1 + r) −r 0 −r 2(1 + r) −r .. .. .. . . . 0 ··· −r 0 ··· 0

 0   0   0  ..  .   2(1 + r) −r  0 1 ··· ··· ··· .. .

Die Bestimmung der Werte f¨ur einen Zeitschritt u j+1 mit Hilfe der vorher berechneten Werte u j ergibt sich dann zu u j+1 = C−1 Au j Das Crank-Nicolson-Verfahren ist f¨ur alle r > 0 stabil und l¨aßt sich leicht auf drei oder mehr Dimensionen erweitern (nat¨urlich kann man es dann graphisch nicht mehr so sch¨on darstellen). Wir wollen hier aber nicht weiter auf Einzelheiten eingehen und verweisen auf die folgende Literatur: • C. Gerald, P. Wheatley, Applied Numerical Analysis, 7. Auflage (2004), Pearson Education, Inc. • A. Hoffmann, B. Marx, W. Vogt, Mathematik f¨ur Ingenieure 2, 1. Auflage (2006), Pearson Studium • G. Venz, L¨osung von Differentialgleichungen mit programmierbaren Taschenrechnern, 2. Auflage (1980!!!), Oldenbourg Verlag

4.2

Finite-Elemente-Methode (FEM)

Ein großer Nachteil der FDM ist, dass sie zur Formulierung der endlichen Differenzenquotienten als Ersatz f¨ur die exakten Differentialquotienten ein kartesisches Gitter (beliebiger Dimension) ben¨otigt. Dies f¨uhrt unweigerlich zu einem Problem, wenn es darum geht komplizierte, z.B gekr¨ummte bzw. detailreiche R¨ander eines Gebiets zu erfassen. In Abb. 11 ist so eine Situation dargestellt. Einen optisch“ a¨ hnlichen, aber konzeptionell v¨ollig anderen Ansatz verfolgt die Finite-Elemente-Methode ” (FEM). Bei der FEM wird das Problemgebiet in kleine, durch Strecken oder B¨ogen berandete Gebiete (z.B. Geraden (1D), Dreiecke und Vierecke (2D), Tetraeder und Quader(3D)) aufgeteilt. Man spricht von einer l¨uckenlosen Pflasterung“ oder Triangulierung des Gebiets mit sog. Finiten Elementen. Die L¨osung des Problems ”

Angewandte Mathematik II

24

Abbildung 11: Gebietsdiskretisierung bei einem a) FDM und einem b) FEM Verfahren.

wird dann st¨uckweise mit Hilfe einfachen Funktionen, die in jedem Element definiert sind, n¨aherungsweise bestimmt, vergleichbar mit der Spline-Interpolation. Die unbekannten Koeffizienten dieser Element-Funktionen werden mit Hilfe einer Variationsformulierung f¨ur die PDE, inklusive ihrer Anfangs- und Randbedingungen bestimmt. Dies f¨uhrt wie im Fall der FDM auf ein (ebenfalls d¨unn besetztes) Gleichungssystem zur Bestimmung der Funktionswerte an den Ecken“ der Finiten Elemente (und auch innerhalb dieser). Gemeinsam mit ” den Element-Funktionen ergibt sich daraus eine st¨uckweise zusammengesetzte L¨osung auf dem gesamten Gebiet. Die FEM ist deutlich komplexer als die FDM, aber wesentlich m¨achtiger. Wir werden die Methode in den Folgenden Abschnitten anhand unseres beliebten Beispiels der W¨armeleitungsgleichung in einer Dimension entwickeln. Der Vorteil hierbei ist, dass jeder Schritt per Hand“ nachvollzogen werden kann und die Metho” de einmal mit Bleistift und Papier“ durchgef¨uhrt wird. Damit ergibt sich ein grunds¨atzliches Verst¨andnis der ” zugrunde liegenden Algorithmen“, die in den unz¨ahligen FEM Werkzeugen implementiert sind. Die eindimen” sionale Variante ist v¨ollig ausreichend f¨ur das Verst¨andnis der FEM; die Erweiterung auf h¨ohere Dimensionen und kompliziertere PDE erscheint dann offensichtlich, und zwar als komplexe, aber nicht komplizierte Aufgabe. 4.2.1

Entwicklung der Finiten Elemente Methode

Bei der FEM handelt es sich um ein modernes Verfahren, das haupts¨achlich in den letzten 50 Jahren entwickelt wurde. Die Vorgeschichte der FEM und ihrer Grundlagen beginnt jedoch viel fr¨uher. Im 19. Jahrhundert beschreibt Schellbach8 die L¨osung eines Minimalfl¨achenproblems, das der FEM schon recht a¨ hnelte. Weitere Meilensteine waren die 1909 bzw. 1915 erschienenen Aufs¨atze von Ritz9 und Galerkin10 . Die Nutzbarmachung der Variationsrechnung kam jedoch erst Mitte des 20. Jahrhunderts, als man erkannte, dass Problemstellungen viel leichter l¨osbar sind, wenn man lokale Ansatzfunktionen z.B. zur Beschreibung der Verbiegung eines Balkens durch ein Variationsprinzip angeht. Bis dahin durchgef¨uhrte Untersuchungen auf Basis von Variationsprinzipien verwendeten globale Ansatzfunktionen, die sich u¨ ber die gesamte Struktur erstreckten. Diese Formulierungen gehen auf Courant11 zur¨uck, der beschreibt, wie im Ritzschen Verfahren Ansatzfunktionen mit lokalem Tr¨ager zu verwenden sind, und wie FEM Berechnungen nach diesem Prinzip durchgef¨uhrt werden. Courants Arbeit konnte jedoch damals keine praktische Bedeutung erlangen, da es bei der Verwendung von elementweise begrenzten Ansatzfunktionen zu großen Gleichungssystemen kommt, die ohne elektronische Rechenanlagen noch nicht l¨osbar waren. Die Notwendigkeit zur Untersuchung von Tragwerken mit 50 und mehr Freiheitsgraden ergab sich zuerst im Flugzeugbau. Die M¨oglichkeiten zum Einsatz geeigneter Rechner bot sich ab ca. 1950. Dazu war es jedoch zuerst notwendig, die Statik durch Einf¨uhrung der Matrixschreibweise klar und u¨ bersichtlich zu gestalten, um 8 Karl

Schellbach: Variationsrechnung; L¨osung eines Minimalfl¨achenproblems (1851/52) Ritz (1878-1909): neue Methode zur L¨osung von Variationsproblemen, Ritzsches Verfahren (1908/09) 10 Boris G. Galerkin: Verfahren der gewichteten Residuen (1915) 11 Richard Courant (1888-1972): Variational methods for the solution of problems of equilibrium and vibration(s). 1943 (Ansatzfunktionen mit lokalem Tr¨ager, elementweise Ans¨atze f¨ur Schwingungsprobleme) 9 Walter

Angewandte Mathematik II

25

die Rechenabl¨aufe zu schematisieren. Ein entscheidender Verdienst kommt hierbei John H. Argyris12 zu, dessen Arbeiten an der Universit¨at Stuttgart den entscheidenden Beitrag zur kommerziellen Entwicklung der FEM darstellen. Die erste, Struktur- und Str¨omungsmechanik umfassende, Anwendung der FEM, die als Geburts” urkunde“ der FEM gilt, wurde 1956 bei Boeing an gepfeilten Flugzeugfl¨ugeln durchgef¨uhrt13 . Die Entwicklung der FEM ist bei weitem noch nicht abgeschlossen. Heutige Trends im Computer Aided Engineering (CAE) sind • die Verkn¨upfung von CAD/CAM mit der FE-Analyse, • die Vereinheitlichung von Netzen f¨ur unterschiedliche Simulationsschritte (Beispiel: Blechumformung, Strukturberechnung, Aeroelastizit¨at bis hin zur Crashsimulation mit einem Netz, Stichwort: Unified ” FEA“), • objektorientierte Modelle, • die Prozessautomatisierung innerhalb der FE-Pakete und • gekoppelte Simulationen (Stichwort: Multiphysics“). ” 4.2.2

Problemformulierung

Am Beginn steht wie in allen anderen Verfahren die Problemformulierung (engl. problem statement)14 . Bevor wir mit dem konkreten Beispiel beginnen, wollen wir zun¨achst einige allgemeine Zusammenh¨ange der Methode erl¨autern. Grunds¨atzlich erfolgt die Problemformulierung wie gehabt als PDE (-System) f¨ur die gesuchte, unbekannte Gr¨oße u = u(~x). Man kann dies abk¨urzend schreiben als L[u] = 0

in Ω ⊂ ’n

(4.1)

Hierbei ist L ein Platzhalter f¨ur die PDE; Die unabh¨angigen Variablen seien zu einem Vektor ~x = (x, y, z, t)T zusammengefasst. Beispiel: L[u] = ut − a2 u xx Die Anfangs- und Randbedingungen seien abk¨urzend dargestellt als l(u) = 0

auf: ∂Ω ⊂ ’

(4.2)

Beispiel: l(u) = l(u) =

u(0, t) = u(L, t) = 0 x  u(x, 0) = sin π L

Die exakte analytische L¨osung u(~x) von (4.1) mit (4.2) ist f¨ur realistische Probleme praktisch nicht berechenbar! Der entscheidende Ansatz der FEM zur Berechnung einer N¨aherungsl¨osung f¨ur u(~x) lautet: u(~x) ≈ uN (~x) ≡ |{z} N¨aherung!

N X

ψα (~x) · Qα

(4.3)

α=1

Welche Idee steckt hinter der Gleichung (4.3) und damit der FEM? Die N¨aherung uN (~x) ist offensichtlich eine endliche Reihenentwicklung bzw. eine Linearkombination von N Funktionen ψα (~x) mit Koeffizienten Qα (Zahlen!). Im Gegensatz zu anderen Verfahren beginnt die FEM mit 12 John

Argyris (1913-2004): Kraft- und Verschiebungsmethode f¨ur Stabtragwerke, Matrizenformulierung (1954/55). J. Turner, Ray W. Clough, H. C. Martin, L. J. Topp: Stiffness and deflection analysis of complex structures (1956). 14 Die Verwendung englischer Termini soll den Bezug zur u ¨ berwiegend englischen Literatur aber insbesondere zu den FEM-Werkzeugen bzw. deren Handb¨uchern erleichtern. 13 M.

Angewandte Mathematik II

26

einer n¨aherungsweisen Darstellung der noch unbekannten L¨osung. Sie gibt quasi ein Korsett“ f¨ur diese N¨ahe” rungsl¨osung vor. Wie sp¨ater noch erl¨autert wird spannen die ψα (~x) einen N-dimensionalen Vektorraum auf, der die L¨osung tr¨agt“. Man nennt die ψα daher auch Ansatzfunktionen (engl. trial functions). Die Ansatzfunktio” nen ψα werden sp¨ater geschickt“ gew¨ahlt, also vorgegeben. Die unbekannten Koeffizienten Qα werden unter ” Verwendung von (4.1) und (4.2) berechnet. Man erh¨alt eine n¨aherungsweise L¨osung des Problems in der Form (4.3)). Wie berechnet man die Koeffizienten Qα , und zwar so, dass die N¨aherungsl¨osung nur einen m¨oglichst kleinen Fehler gegen¨uber der exakten L¨osung hat? Die entscheidende Idee hierzu liefert die Variationsrechnung, speziell die Methode der gewichteten Residuen. 4.2.3

Variationsformulierung - Die Methode der gewichteten Residuen

Setzt man die exakte L¨osung u(~x) in die PDE (4.1) ein, so gilt selbstverst¨andlich: L[u] = 0 Angenommen man h¨atte bereits auf irgendeine Weise die N¨aherungsl¨osung, wie sie in Gleichung (4.3) definiert ist, bestimmt. Was passiert, wenn man diese N¨aherungsl¨osung uN (~x) in die PDE (4.1) einsetzt? Nun, da die L¨osung nicht exakt ist, ergibt sich ein Fehler bzw. es bleibt ein von Null verschiedener Rest u¨ brig, also h i h i L uN , 0 bzw. L uN = R(~x) , 0. Man nennt R(~x) das Residuum (lat. f¨ur Rest) und er quantifiziert den Fehler der N¨aherung. Das Ziel ist nun, die Qα so zu bestimmen, dass das Residuum R(~x) m¨oglichst klein wird. Hat man das erreicht, kann man guter Hoffnung sein, dass die N¨aherung nicht allzu weit von der wahren L¨osung entfernt ist (auch diese Aussage kann man quantifizieren!). F¨ur solch eine Optimierungsaufgabe gibt es in der Mathematik ein Fachgebiet, dass sich Variationsrechnung nennt. Innerhalb der Variationsrechnung gibt es ein Verfahren, dass sich Methode der ” gewichteten Residuen“ nennt und hier angewendet werden soll bzw. muss. Wie geht man bei diesem Verfahren vor? Man multipliziert zun¨achst das Residuum R(~x) mit sogenannten Gewichts- oder Testfunktionen (engl. test functions) ϕβ (~x) und fordert: Z ! ϕβ (~x) R(~x) d~x = 0 ∀ ϕβ (~x), 1 ≤ β ≤ N (4.4) Ω

Wichtig ist hier der Zusatz ; ∀ ϕβ (~x), d.h. wir haben hier nicht eine Gleichung, sondern N Gleichungen f¨ur jede der N Testfunktionen ϕβ (~x) eine. Warum? Nun, im Ansatz f¨ur die N¨aherungsl¨osung Gl. (4.3) gibt es N unbekannte Koeffizienten Qα . Um diese eindeutig zu bestimmen, ben¨otigt man genau N Bedingungen! Diese liefern eben die Gleichungen (4.4). Durch die Integration verschwindet die ~x-Abh¨angigkeit und es bleibt ein einfaches algebraisches Gleichungssystem von N Gleichungen zur Bestimmung der N unbekannten Qα u¨ brig! Warum nun Gleichung (4.4) ausgerechnet diejenigen Werte f¨ur die Qα liefert, die den Fehler der N¨aherung Gl. (4.3) minimieren, soll sp¨ater erl¨autert werden. Wir wollen uns zun¨achst der Frage widmen wie die Ansatzund Testfunktionen ψα und ϕβ aussehen m¨ussen. Dazu verwenden wir ein einfaches begleitendes Beispiel. Der Ausgangspunkt ist wieder einmal die W¨armeleitungsgleichung. Wir wollen zus¨atzlich eine W¨armequelle s (Einheit W/m3 ) im Gebiet Ω ber¨ucksichtigen. Wir k¨onnten damit z.B. einen Heizstab oder einen Tauchsieder modellieren, der gleichm¨aßig und homogen u¨ ber das Gebiet Ω W¨arme an die Umgebung abgibt. In Comsol oder anderen FEM-Werkzeugen findet man zur Modellierung der W¨armeleitung folgende Formulierung der W¨armeleitungsgleichung:   ~ k∇u(~ ~ x, t) + ∂u(~x, t) = s in Ω ⊂ ’3 −∇ ∂t | {z } ?

Angewandte Mathematik II

27

Man kann den Term ? als Rate, mit der Energie durch ein Medium mit ortsabh¨angiger Temperaturleitf¨ahigkeit k transportiert wird, interpretieren. Beachten Sie, dass die Ortsabh¨angigkeit von k dazu f¨uhrt, dass man k nicht nach vorne ziehen und die beiden Nabla-Operatoren zu einem Laplace-Operator zusammenfassen kann. Erst, wenn k als konstant angenommen wird, darf man das tun und erh¨alt z.B. Gl. (4.1) aus dem vorherigen Abschnitt. Wie in den vorherigen Beispielen wollen wir uns auch hier auf den eindimensionalen und station¨aren (=zeitunabh¨angigen) Fall beschr¨anken: u = u(x) ;





! d du k −s = dx dx ⇔ L[u] =

∂u = 0; ∂t

~ → ∂ → d { ∇ ∂x dx

0 0

(4.5) (vgl. Gl.(4.1))

im Gebiet Ω : 0 ≤ x ≤ L

Zur L¨osung dieser station¨aren W¨armeleitungsgleichung ben¨otigen wir keine Anfangsbedingungen, wohl aber Randbedingungen. Wir verwenden: An der Stelle x = 0: W¨armefluß f n : −k

du = fn dx

(4.6)

({ Neumann Randbedingung), An der Stelle x = L: Konstante Temperatur:

u(x = L) = T b

(4.7)

({ Dirichlet Randbedingung) F¨ur einen sp¨ateren Vergleich geben wir die exakte L¨osung des Problems an: "  x 2 # f n L  s L2 x u(x) = 1− 1− + Tb + 2k L k L

(4.8)

Gem¨aß Gl. (4.3) verwenden wir f¨ur die N¨aherungsl¨osung den Ansatz u(x) ≈ uN (x) =

N X

ψα (x) Qα ,

(4.9)

α=1

wobei die Ansatzfunktionen ψα (x) entsprechend der Problemformulierung nun eindimensional sind. F¨ur die Variationsformulierung ben¨otigen wir das Residuum, das in unserem Fall konkret lautet: ! h i d duN N R(x) = L u (x) = − k − s. (4.10) dx dx Setzt man diesen Ausdruck in Gleichung (4.4) ein, so ergibt sich Z ! ϕβ (x) R(x) dx = 0

∀ ϕβ (x)

bzw. 1 ≤ β ≤ N



Z ⇔ Ω

! # d duN (x) ϕβ (x) − k − s dx dx dx "

!

=

0;

1 ≤ β ≤ N.

(4.11)

Angewandte Mathematik II

28

Dies sind die N Gleichungen zur Bestimmung der N unbekannten Koeffizienten Qα . Als n¨achstes setzen wir den Ansatz uN (x) ein und erhalten    N   Z    d  d X 1 ≤ β ≤ N. (4.12) ϕβ (x) − k  ψα (x) Qα  − s dx = 0 ; dx dx α=1 Ω

Da die Integration eine lineare Operation ist, d¨urfen wir Summation und Integration vertauschen:   !# " Z Z   N   X   d dψ (x)   α ϕ (x) − Qα ·  k dx − ϕ (x) s dx = 0; 1 ≤ β ≤ N.  β β     dx dx   α=1 Ω

(4.13)



Dies ist der erste entscheidende Schritt, da nun die zu bestimmenden Koeffizienten Qα aus den Integralen herausgezogen worden sind. Beim Blick auf Gl. (4.13) f¨allt auf, dass die Ansatzfunktionen ψα (x) zweimal differenzierbar sein m¨ussen, w¨ahrend die Testfunktionen ϕβ (x) gar nicht differenzierbar sein m¨ussen! Dies ist eine starke“ Forderung an ” die Wahl der ψα (x)! Warum? Man schließt mit dieser Forderung sofort lineare Ansatzfunktionen aus. Interessanterweise kann man durch eine einfache Umformung diese Forderung abschw¨achen“ und zwar durch ” Anwendung der partiellen Integration! Zb Zur Erinnerung:

 b Zb f (x) g (x) dx = − f (x) g(x) + f 0 (x) g(x) dx 0



a

a

a

 dψβ (x)  d k dx f¨uhrt die Anwendung der partiellen Integration auf Mit f (x) = ψβ (x) und g0 (x) = − dx   ! ! x=L  Z Z N  X  dϕ (x) dψ (x) dψ (x) β α α Qα ·  − ϕβ (x) s dx = 0 ; k dx − ϕβ (x) k 1 ≤ β ≤ N,  x=0  dx dx dx α=1

(4.14)





und wir erhalten damit unsere gew¨unschte Abschw¨achung der Forderung nach Differenzierbarkeit. Die Ansatzfunktion ψα (x) m¨ussen nur noch einmal differenzierbar sein, daf¨ur allerdings die Testfunktionen aber auch. Gl.(4.14) wird deshalb auch schwache Formulierung“ (engl. weak statement) genannt! ” Obwohl wir die schwache Formulierung“ an einem konkreten Beispiel entwickelt haben, so ist doch die Vor” gehensweise v¨ollig Problem- bzw. PDE-unabh¨angig. Andererseits gibt es die schwache Formulierung“ aber ” nicht als universelle Formel f¨ur alle Problemformulierungen, sie muss f¨ur jede PDE individuell hergeleitet werden. Die schwache Formulierung“ ist jedoch v¨ollig unabh¨angig von der Wahl der Ansatz- und Testfunk” tionen ψα und ϕβ . Und es ist in der Tat so, dass sich je nach Definition dieser Ansatz- und Testfunktionen die unterschiedlichen numerische Verfahren ergeben, die auf dem Markt erh¨altlich sind, wie z.B. die FEM, aber auch die BEM (Boundary Element Method) und andere weniger verbreitete Verfahren. In Vorausschau auf die Entwicklung der FEM in den kommenden Abschnitten ist noch eine weitere Eigenschaft der schwachen For” mulierung“ Gl.(4.14) wichtig. Sie gilt n¨amlich auch f¨ur nur st¨uckweise glatte Funktionen. Auch, wenn die schwache Formulierung“ zun¨achst unabh¨angig von den Ansatz- und Testfunktionen ent” wickelt wurde, so stellt sich die Frage, ob es nicht optimale Ansatz- und Testfunktionen gibt? Was heißt in diesem Zusammenhang optimal“? Optimal sind diejenigen Ansatz- und Testfunktionen, die den Fehler der ” N¨aherung uN (x) minimieren! Interessanterweise kann man den Fehler mit Hilfe einer geeigneten Norm (darauf kommen wir sp¨ater noch) genau dann minimieren, wenn die Ansatz- und Testfunktion identisch sind, also gilt: ψα = ϕ β ;

α=β;

1≤α≤N

(4.15)

Der Beweis hierf¨ur wurde von Galerkin erbracht und es sei darauf hingewiesen, dass es sich bei diesem Trick“ ” nicht um eine N¨aherung handelt!

Angewandte Mathematik II

29

4.2.3.1 Galerkin-Formulierung Mit Gleichung (4.15) kann man die schwache Formulierung“ Gl. (4.14) umgeschrieben zu: ”   ! ! x=L  Z Z N  X  dψ (x) dψ (x) dψ (x) β α α Qα ·  − ψβ (x) s dx = 0 ; 1≤β≤N k dx − ψβ (x) k  x=0  dx dx dx α=1

(4.16)





Man nennt Gl. (4.16) die Galerkin Formulierung“ (engl. Galerkin weak statement (GWS)). Was beim Betrach” ten von Gl. (4.16) sofort ins Auge f¨allt, ist dass Fluß-Randbedingungen quasi automatisch“ eingebaut sind! Es ” gilt n¨amlich   ! ! x=L N  X   duN dψα (x)    n − ψ (x) k = −ψ (x) f −ψ (x) k (4.17) Qα   β β     β dx dx x=L x=0 x=0 α=1

Die Galerkin-Formulierung erleichtert uns nun nochmals das Leben, da wir jetzt nur noch einen Satz von Funktionen {ψα } finden m¨ussen. Wie sehen die {ψα } nun aus? Noch einmal zur Erinnerung: Sobald die {ψα } bekannt sind, k¨onnen die Qα berechnet werden und das Problem ist (n¨aherungsweise) gel¨ost! Die Mindestanforderung an die Ansatz- und Testfunktionen {ψα } kann man der Galerkin-Formulierung Gl. (4.16) direkt entnehmen. Die dort auftauchenden Integrale m¨ussen zun¨achst einmal u¨ berhaupt existieren, d.h. es muß gelten: Z dψ (x) dψ (x) β α dx < ∞ (4.18) dx dx Ω

Diese Forderung ist mathematisch a¨ quivalent zu der Aussage, dass die {ψα } den Hilbert-Raum H 1 aller Funktionen aufspannen, deren erste Ableitung quadratintegrabel ist! 4.2.4

Gebietsdiskretisierung - Finite-Elemente-Ansatzfunktionen

Die Galerkin-Formulierung Gl. (4.16) hat bislang nichts mit Finiten Elementen zu tun! Die Ansatz- bzw. Testfunktionen {ψα (x)} sind immer noch global auf dem gesamten Gebiet Ω (hier, der Stab mit Ω : 0 ≤ x ≤ L) definiert. Die Idee hinter der FEM ist ja aber gerade die, dass man anstatt nach einer L¨osung, die sich als Linearkombination globaler Funktionen {ψα (x)} darstellen l¨aßt, versucht diese st¨uckweise mit Hilfe lokaler Funktionen zu berechnen, die dann aneinander geheftet werden. Das Problem der Darstellung einer Funktion mit Hilfe globaler Ansatzfunktionen kennt man beispielsweise aus der Fourier-Synthese oder der Polynominterpolation. Um eine m¨oglichst gute N¨aherung der gesuchten Funktion mit globalen Ansatzfunktionen im gesamten Gebiet Ω zu bekommen, ben¨otigt man in der Regel sehr viele (N → ∞) und komplexe Ansatzfunktionen (eben z.B. Sinus- und Cosinus-Funktionen), um alle lokalen Feinheiten erfassen zu k¨onnen. In manchen F¨allen reicht aber selbst das nicht aus, um pathologische“ ” ¨ Funktionsverl¨aufe ohne Probleme darstellen zu k¨onnen15 Mit Hilfe stuckweise definierter lokaler Ansatzfunktionen kann man eine gute N¨aherung auch mit nur wenigen (N  ∞) und dar¨uber hinaus auch einfachen Ansatzfunktionen (z.B. Polynome niedriger Ordnung) realisieren. Man kann diesen Paradigmenwandel mit dem ¨ Ubergang von der (globalen) Polynominterpolation zur (lokalen) Spline-Interpolation vergleichen. In Abb. 12 ist das Problem skizziert. Die Darstellung einer bekannten Funktion f (x) mit Hilfe eines globalen Interpola¨ Die Darstellung mit Hilfe lokaler, st¨uckweiser, kubischer tionspolynoms f¨uhrt zu heftigen Uberschwingern“. ” Splines hingegen entspricht fast perfekt der vorgegebenen Funktion. Warum ist das so? Die Ursache liegt im Fundamentalsatz der Algebra (siehe Angewandte Mathematik I“). Danach muss ein Polynom seinem Grad ” entsprechend viele Nullstellen besitzen. Der Grad des Interpolationspolynoms ergibt sich aus der Anzahl der 15 Man denke nur an die unvermeidlichen Uberschwinger“ ¨ bei der Darstellung eines Rechteckpulses mit Hilfe von Sinus- und Cosinus” Funktionen, dem sog. Gibbschen Ph¨anomen.

Angewandte Mathematik II

30

Interpolationspunkte bzw. St¨utzstellen. Beim globalen Interpolationspolynom m¨ussen alle St¨utzstellen ber¨ucksichtigt werden, d.h. der Grad des Polynoms ist entsprechend hoch. F¨ur die st¨uckweise definierten Splines wird aber immer nur zwischen zwei St¨utzstellen interpoliert. Hier reichen grunds¨atzlich Polynome 3. Grades aus, um bestm¨oglich (mit der geringsten Kr¨ummungsenergie“) zu interpolieren, und zwar unabh¨angig von der Anzahl ” globaler St¨utzstellen! Wir halten fest, dass die st¨uckweise Darstellung einer Funktion folgendes erfordert:

Abbildung 12: Darstellung einer bekannten Funktion f (x) mit Hilfe eines globalen Interpolationspolynoms und mit Hilfe von Spline-Funktionen (Polynome 3. Ordnung).

I Eine Aufteilung des globalen Gebiets in lokale Gebiete ( Elemente“), in denen nur wenige St¨utzstellen ” (auch Knoten genannt, engl. nodes) liegen. II Eine einfache (Interpolations-) Funktion je lokalem Gebiet, deren Verlauf nur durch die lokalen St¨utzstellen festgelegt ist sowie III Forderungen an die Glattheit und Differenzierbarkeit der Vereinigungsmenge aller lokalen Funktionen, so dass die globale Darstellung mindestens stetig verl¨auft. Gem¨aß dieser Punkte ist die einfachste aller Interpolationen diejenige mit Hilfe von Geradest¨ucken wie in Abb.13 gezeigt! Beachten Sie, dass der Verlauf zwar u¨ berall stetig, aber nicht u¨ berall differenzierbar ist! 4.2.4.1 Veranschaulichung via 2-Elemente-Gebietsdiskretisierung Wir wollen nun die Finite-Elemente-Formulierung f¨ur unsere Problemformulierung mit Hilfe der eindimensionalen W¨armeleitungsgleichung an einem einfachen Beispiel, n¨amlich einer Gebietsdiskretisierung mit zwei Elementen veranschaulichen (siehe Abb. 14). Wir nennen die Elemente Ωe=1 und Ωe=2 . An den Enden der Elemente befinden sich die sog. Knoten mit der Bezeichnung X1, X2, X3. Bei der Betrachtung der Gebietsdiskretisierung stellt sich eine fundamentale Frage: Wie l¨aßt sich der Ansatz PN f¨ur die N¨aherungsl¨osung uN (x) = α=1 ψα (x) Qα mit der Element-weisen Gebietsdiskretisierung verkn¨upfen? Das fundamentale an dieser Frage ist, dass ihre Beantwortung unmittelbar eine geometrische Interpretation des Ansatzes liefert! N¨amlich die Folgende: • Die Koeffizienten Qα werden als Funktionswerte an den Knoten Xα interpretiert. • Die Ansatzfunktionen ψα dienen der Interpolation zwischen den Knoten.

Angewandte Mathematik II

31

f(x) u(x) (exakt) N

u (x)

x1

0

x2

x3

L

x

Abbildung 13: Die einfachste, nicht aber unbedingt die beste, Interpolation zwischen St¨utzstellen wird durch Geradest¨ucke (Polynome 1. Grades) realisiert. Q Q

001 11 11 00 00 11

Ω e=1

Q

2 11 00 00 11 00 11

Ωe=2

X2

X1 (0)

Q

3 11 00 00 11 00 11

X3 (L)

Abbildung 14: Veranschaulichung der FEM mit Hilfe einer 2-Elemente-Gebietsdiskretisierung zur n¨aherungsweisen L¨osung des W¨armeleitungsproblems in einer Dimension.

In unserem konkreten Beispiel bedeutet das f¨ur die N¨aherungsl¨osung u (x) = N

N X

ψα (x) Qα = ψ1 (x) Q1 + ψ2 (x) Q2 + ψ3 (x) Q3 .

α=1

Eine weitere wichtige Eigenschaft der N¨aherungsl¨osung wird mit Hilfe der geometrischen Interpretation deutlich: Nur die Qα haben eine physikalische Einheit (hier z.B. Kelvin); die Ansatzfunktionen sind einheitslos. Beachten Sie, dass die Ansatzfunktionen ψα (x) immer noch global definiert sind, also Ω : 0 ≤ x ≤ L gilt. Im ersten Schritt wollen wir das auch beibehalten und uns mit Hilfe der Anforderungen an die N¨aherungsl¨osung die konkrete Form der Ansatzfunktionen beschaffen. Gem¨aß geometrischer Interpretation muß ... I uN (x) an der Stelle X1 den Wert Q1 haben, d.h. uN (X1) = ψ1 (X1) Q1 + ψ2 (X1) Q2 + ψ3 (X1) Q3 , woraus folgt ψ1 (X1)

= 1

ψ2 (X1) = ψ3 (X1)

= 0,

II uN (x) an der Stelle X2 den Wert Q2 haben, d.h. uN (X2) = ψ1 (X2) Q1 + ψ2 (X2) Q2 + ψ3 (X2) Q3 ,

Angewandte Mathematik II

32

woraus folgt ⇒ ψ1 (X2) = ψ3 (X2)

= 0

ψ2 (X2)

= 1,

III uN (x) an der Stelle X3 den Wert Q3 haben, d.h. uN (X3) = ψ1 (X3) Q1 + ψ2 (X3) Q2 + ψ3 (X3) Q3 , woraus folgt ⇒ ψ1 (X3) = ψ2 (X3)

= 0

ψ3 (X3)

= 1.

Q Q

1 11 00 00 11 00 11

X1 (0)

ψ1

Ω e=1

Q

2 11 00 00 11 00 11

Ωe=2

X2

Q3 1 0 0 1

X3 (L)

11 00 00 11 00 11

11 00 00 11 00 11

ψ2

11 00 00 11 00 11

11 00 00 11 00 11

11 00 00 11 00 11

00 11 11 00 00 11

11 00 00 11 00 11

ψ3 00 11 11 00 00 11

00 11 11 00 00 11

Abbildung 15: Die Ansatzfunktionen ψi (x), hergeleitet mittels geometrischer Interpretation. Damit kennen wir die Werte aller drei Ansatzfunktionen ψα (x) an den Stellen Xα! Was fehlt sind die Werte der ψα (x) zwischen den Knoten, d.h. die Interpolation. Im einfachsten Fall verbinden wir einfach die Knotenwerte mit Geradenst¨ucken, sprich Polynomen 1. Grades. Damit erhalten wir die vollst¨andige Kenntnis aller Ansatzfuntkionen ψα (x) wie in Abb. 15) dargestellt. Wie w¨urde sich das Bild a¨ ndern, wenn wir statt zwei, drei Elemente und damit vier Ansatzfunktionen verwenden? Die Systematik ist unmittelbar erkennbar. Die ψα (x) sehen aus wie ein wanderndes Dreieck“ 16 , das sich ” u¨ ber drei Knoten erstreckt, wobei die Spitze von ψα an der Stelle Xα ist und die Schenkel von X(α − 1) zu Xα sowie von Xα zu X(α + 1) verlaufen! Dieses Ergebnis ist v¨ollig unabh¨angig von der Anzahl der Elemente bzw. der Anzahl der Knoten N, gilt also beispielsweise auch f¨ur 150 Knoten! Das Bild vom wandernden Dreieck“ ” 16 In

der FEM Literatur wird diese Ansatzfunktion auch Hut-Funktion“ genannt. ”

Angewandte Mathematik II

33

l¨aßt sich auch mathematisch formulieren:  x−X(α−1)   , f¨ur X(α − 1) ≤ x ≤ Xα  X(α)−(α−1)    X(α+1)−x ψα (x) =  , f¨ur Xα ≤ x ≤ X(α + 1)  X(α+1)−Xα    0 , sonst, d.h. f¨ur x < X(α − 1) ; x > X(α + 1)

(4.19)

Hierbei sind zwei Dinge auff¨allig: I Jede Ansatzfunktion ψα ist nur in der unmittelbaren Umgebung der Stelle Xα von Null verschieden, also stark lokalisiert. II Es handelt sich trotzdem um globale“ Funktionen, die auf dem gesamten Gebiet 0 ≤ x ≤ L definiert ” sind wie die letzte Zeile verr¨at. Noch einmal zur Erinnerung. Die so definierten Ansatzfunktionen tragen“ die N¨aherungsl¨osung uN (x) und ” ¨ zwar in der Summe bzw. als Uberlagerung. Stellen wir deswegen alle drei ψα (x) einmal in einem Graphen dar, so ergibt sich eine Darstellung wie in Abb. 16 gezeigt. Was f¨allt in dieser Darstellung auf? Die Projektion

u2

u1

1

11 00 00 11 00 11

X1

Ω e=1

11 00 00 11 00 11

X2

Ωe=2

11 00 00 11 00 11

X3

Abbildung 16: Darstellung aller Ansatzfunktionen in einem Graphen; lokale Funktionen u1 und u2 in einen Graphen hat zu einem J¨agerzaunmuster“ im gesamten Gebiet gef¨uhrt. Das mag vielleicht a¨ sthetisch ” umstritten sein, aber es zeigt eines deutlich: Es gibt keinen Unterschied zwischen der Darstellung in Element 1 und 2! Und auch hier w¨urde sich am Muster nichts a¨ ndern, wenn man statt zwei Elementen 150 verwendet! Wir k¨onnen jetzt unmittelbar zu einer lokalen Darstellung u¨ bergehen. Angenommen wir wollen das J¨agerzaunmu” ster“ mit Hilfe lokaler Funktionen, die in jedem Element definiert sind, erzeugen. Wie viele Funktionen ben¨otigt man dazu? Tats¨achlich nur zwei und nicht mehr! Diese zwei Funktionen u1 und u2 sind in beiden Elementen identisch! Was w¨urde sich a¨ ndern, wenn man statt zwei Elementen 150 verwendet? Nichts! D.h. auch in dem Fall ben¨otigt man nur zwei lokale Elementfunktionen u1 und u2 ! In einer globalen Definition wie oben w¨urde man allerdings 151 (Anzahl der Knoten=Anzahl der Elemente plus eins!) Ansatzfunktionen ψα (x)ben¨otigen! ¨ Der Ubergang zu einer st¨uckweisen, lokalen Darstellung f¨uhrt somit zu einer großen Vereinfachung des Problems – genau wie im oben diskutierten Fall der Polynominterpolation. 4.2.4.2 Finite-Elemente-Ansatzfunktionen Da die FE-Ansatzfunktionen in jedem Element identisch sind, kann man sie allgemein angeben. Es gilt u1 (x) =

XR − x XR − XL

und

u2 (x) =

x − XL XR − XL

Hierbei sind • XR und XL die lokalen R¨ander“ irgendeines Elements, ” • XR − XL = le ist die elementspezifische Breite.

in Ωe : XL ≤ x ≤ XR

(4.20)

Angewandte Mathematik II

34

Beachten Sie, in dieser lokalen Darstellung werden lokale (Element-) Koordinaten verwendet. Wir m¨ussen als an irgendeiner Stelle die Verbindung zwischen dieser Darstellung und dem globalem L¨osungsgebiet 0 ≤ x ≤ L herstellen! F¨ur das weitere Vorgehen, verwenden wird die folgende kompakte Darstellung der FEAnsatzfunktionen als Vektor: ( ) 1 XR − x (4.21) {N} = le x − X L e Im Gegensatz zu den globalen Ansatzfunktionen ψα (x) aus Gleichung (4.19) sind die FE-Ansatzfunktionen (4.21) nun lokal je Element Ωe definiert. Damit k¨onnen wir die globalen Integrale der Galerkin-Formulierung (4.16) u¨ ber das Gesamtgebiet Ω nun in eine Summe u¨ ber lokale Integrale u¨ ber die Elemente Ωe des FiniteElement-Gitters ersetzen. 4.2.5

Finite-Elemente-Formulierung

Um zu einer FE-Formulierung der Galerkin-Formulierung zu kommen, m¨ussen wir zun¨achst die N¨aherung mit Hilfe der FE-Ansatzfunktion in eine lokale Elemente-Form u¨ berf¨uhren. Hierzu werden einfach die globalen Ansatzfunktionen durch die lokalen FE-Funktionen ersetzt: ! N M X 2 M X X X Qe1 (u1 u2 ) · uN (x) = (x) ψα (x)Qα = u j (x) · Qej = Qe2 α

=

e=1 j=1

M X

e=1

{N}T · {Q}e

(4.22)

e=1

↑ kein Skalarprodukt! → Matrizenmultiplikation!

Im letzten Schritt wurde eine weitere Vereinfachung vorgenommen, die die Summe u¨ ber die FE-Ansatzfunktion durch die Darstellung als Matrizenmultiplikation ersetzt. Dies f¨uhrt im Folgenden zu einer Vereinfachung der Schreibweise. Damit l¨aßt sich die globale Galerkin-Formulierung N X

Z Qα

α=1



! ! L Z dψβ (x) dψα (x) du N k dx − ψβ · s · dx − ψβ k · =0 dx dx dx 0 Ω

umschreiben zu einer lokalen Galerkin-Formulierung M X 2 X

Z Qej

e=1 j=1

Ωe

! ! x Z dui du j du N R k dx − ui · s · dx − ui k · = 0, dx dx dx xL

1 ≤ i ≤ 2.

(4.23)

Ωe

Die FE-Ansatzfunktionen ui u¨ bernehmen nun die Rolle der (ehemaligen) Ansatz- und Testfunktionen ψβ (x). N Gleichungen hat man nun nur noch zwei Statt x Gleichungen (4.23). Die Auswertung an den globalen R¨andern L wird durch die elementspezifischen R¨ander L ersetzt! 0 P xR Die Summe u¨ ber die FE-Ansatzfunktionen 2j=1 verschwindet bei Verwendung der Matrizenmultiplikation und man erh¨alt die kompaktere Form: 2 X j=1

Z Qej Ωe

! ! Z dui du j d {N} d {N}T k dx = k dx · {Q}e , dx dx dx dx Ωe

( mit: {Q}e =

Qe1 Qe2

) (4.24)

Angewandte Mathematik II

35

Wir wollen nun die einzelnen Terme in Gl. (4.23) f¨ur unser begleitendes Beispiel nacheinander berechnen. Zun¨achst bestimmen wir die folgende Ableitung: ( )! ( ) d {N} d 1 XR − x 1 −1 = = . (4.25) 1 dx dx le x − XL e le ↑ (Gleichung (4.21))

Damit kann der erste Integralausdruck vereinfacht werden zu ! ( ) Z Z 1 −1 1n d {N} d {N}T −1 k dx · {Q}e = ·k· 1 dx le Ωe le Ωe dx !Z 1 1 −1 = kdx · {Q}e . 1 le 2 −1 Ωe

1

o

dx · {Q}e (4.26)

Da die FE-Ansatzfunktionen Geraden sind, und damit ihre Ableitungen konstant, bleibt von der Integration nur noch das Integral u¨ ber k u¨ brig. Nimmt man nun an, dass k konstant sei, so gilt: XR Z Z k dx = k dx = k · x = k · (XR − XL ) = k · le X Ωe Ωe L

und damit 1 le 2

1 −1

−1 1

!Z Ωe

kdx · {Q}e =

k le

1 −1

−1 1

! {Q}e .

Nehmen wir weiterhin an, dass der Quellterm s ebenfalls konstant sei, so ergibt sich: !) Z Z Z ( 1 XR − x dx ui · s · dx , 1 ≤ i ≤ 2 ⇔ {N} · s · dx = s · x − XL Ωe Ωe Ωe le  "( 2 2 2 )# XR  1 s  x · XR − x2  XR 2 − X2R − XR XL + X2L = s· · = x2  XR 2 − X X − X L 2 + X 2 le le  L R L 2 − x · XL XL 2 2 ( 1 ( 1 2 ) ( ) ) 2 s · le 1 s s (XR − XL ) l 2 · 21 e 2 = = = 2 1 1 (X ) l − X le l 2 e R L e 2 2

      

Die Terme, die die Randbedingungen beschreiben, enthalten keine Integrale; hier muss man nur Einsetzen: X X ! X ( ) duN R duN R 1 XR − x duN R ui k · ⇔ k· {N} · k · ; 1≤i≤2 = dx XL dx XL le x − XL dx XL ( ) ( ) duN 1 duN le 1 0 k· − k· = le 0 le dx le dx ( ) ( ) N N du du 0 1 = k· −k· 1 0 dx dx ( ) N du −1 = k· 1 dx Alles eingesetzt in Gl. (4.23) ergibt f¨ur die lokale Galerkin-Formulierung M " X k l=1

le

1 −1

−1 1

!

s · le · {Q}e − 2

1 1

!

duN −k dx

−1 1

!# =0

(4.27)

Angewandte Mathematik II

36

Gl. 4.27 enth¨alt Ausdr¨ucke, die mit dem Index e“ gekennzeichnet sind. Diese Ausdr¨ucke sind elementspezi” fisch und n damitolokal! Beachten Sie, dass die Nummerierung der Knoten in jedem Element gleich ist und zwar immer Qe1 , Qe2 ist. Aber dar¨uber hinaus gilt Gl. (4.27) f¨ur beliebig viele Elemente (der Wert von M ist nicht weiter spezifiziert!)! Wir wollen an unserem begleitendem Beispiel (d.h. mit M = 2) die FE-L¨osung konkret berechnen (siehe L 2

0

11 00 00 11 00 11

Ω e=1

Q1

L

00 11 11 00 00 11

Ωe=2

Q2

00 11 11 00 00 11

Q3

Abbildung 17: Einfachster Fall f¨ur die FE-L¨osung des Beispiels: Aufteilung des Gebiets in zwei Elemente Ω1 und Ω2 mit den globalen Knoten Q1 , Q2 und Q3 . Abb. 17). Das bedeutet, dass wir die L¨osung f¨ur die globalen Knoten Q1 , Q2 , Q3 bestimmen m¨ussen. ABER: Gleichung (4.27) verwendet eine lokale Knoten-Nummerierung Qe1 , Qe2 ! Wie kann man nun die globale und die lokale Nummerierung mit einander in Einklang bringen? Hierzu verwendet man den folgenden Algorithmus“. ” 4.2.5.1 Assemblierung des Finite Elemente Gleichungssystems Ziel ist es, die lokalen {Q}e durch die globalen {Q} zu ersetzen, also etwa ( ⇒

{Q}e=1 =

Q1 Q2

)

( ;

Q2 Q3

{Q}e=2 =

) wird ersetzt durch:

   Q      1   Q2  {Q} =      Q   3

Beim Vergleich f¨allt sofort folgendes Problem auf. In Gleichung (4.27) stehen (korrekterweise!) 2×2 - Matrizen (und 2-Vektoren) f¨ur die elementweise Berechnung der zwei Knoten. Zur Berechnung der 3 globalen Knoten ben¨otigt man aber 3 × 3-Matrizen (und 3-Vektoren)! Die L¨osung des Problems erreicht man durch Erweiterung der 2 × 2 - Matrizen zu 3 × 3 - Matrizen. Hierzu f¨ullt“ man die Spalten und Zeilen der Matrizen (und Vektoren) so mit Nullen auf, dass die globalen Knoten, ” die nicht zum lokalen Element beitragen, auch nicht an der Berechnung beteiligt werden. Diese Operation ergibt f¨ur Element e = 1: • Erster Term: 1 −1

−1 1

!

Q1 Q2

! →

  1   −1 0

   −1 0   Q1     1 0   Q2     0 0 Q3

• Quellterm: 1 1

!

−1 1

!



   1     1  0



   −1     1  0

• Randbedingungen:

und f¨ur Element e = 2:

Angewandte Mathematik II

37

• Erster Term: 1 −1

−1 1

!

Q2 Q3

! →

  0   0 0

  0 0   Q1   1 −1   Q2   Q3 −1 1

   

• Quellterm: 1 1

!

−1 1

!



   1     1  0



   0     −1  1

• Randbedingungen:

Die Summe u¨ ber die Elemente (unter Verwendung der globalen Nummerierung) kann nun durchgef¨uhrt werden. Dieser Schritt wird als Assemblierung (Zusammensetzen; engl. assembly) bezeichnet und f¨uhrt zum folgenden globalen Gleichungssystem          1 −1 0   Q1   1  −1  N    k  du s l      1   1  1 0   Q2  − (4.27) →  −1  1  − k      l1  0 2 dx 0 0 Q3 0 0          0   0 0  0 0   Q1  N    k  du s l       2   −1  = 0 1 −1   Q2  − + (4.28)  1  − k  0    l2 0 −1 2 dx  1  1 1 Q 3

Da beide Elemente die gleiche L¨ange haben gilt:     1 −1 0   2 k     −1 2 −1   ⇒   L  0 −1 1

l1 = l2 = L2  Q1  sL  Q2  −  4 Q3

Man kann Gl. (4.29) umformen zu:     1 −1 0   Q1    2 −1   Q2  −1   0 −1 1 Q3 | {z }

  s L2   = 8k |

A

   1  duN    2  − k dx 1

   −1     0  = 0 1

     1  −1  N    L du    0 ,  2  +   2 dx 1 1 {z }

(4.29)

(4.30)

~b

dem linearen globalen FE-Gleichungsssystem zur Bestimmung der unbekannten Knotenwerte Q1 , Q2 , Q3 . In der FEM-Gemeinde nennt man A die Steifigkeitsmatrix und ~b den Lastvektor. Dies ist eine historisch bedingte Bezeichnung und stammt aus den Anf¨angen der FEM, in denen Verformungsprobleme gel¨ost werden sollten (im Rahmen der lineare Elastizit¨atstheorie). Interessant ist ein Blick auf den letzten Term in (4.29), der die Flußrandbedingungen enth¨alt. Im Rahmen der Umformulierung von der globalen zur elementweisen lokalen Darstellung wurden die Randwerte ebenfalls von global (sprich den Werten 0 und L) zu lokal (XL und XR je Element) umgeschrieben. Nach der Assemblierung bleiben aber tats¨achlich nur die Werte an den globalen Enden (Q1 , Q3 ) von Null verschieden, so wie es auch sein sollte! An den internen Grenzen zwischen den Elementen heben sich die Fl¨usse auf! Einsetzen der Randbedingungen in (4.30): Bei x = 0: −k ·

duN = fn dx

(siehe (4.6)!)

Angewandte Mathematik II Bei x = L:

38

konstante Temperatur T b (siehe (4.7)!)     1 −1 0   Q1    −1 2 −1   Q2    Tb 0 −1 1

  s L2   = 8k

   1  L    2  + 2 k 1

   

fn 0 −F3

    ,

(4.31)

wobei wir abk¨urzend F3 f¨ur den noch unbekannten Fluss am Knoten Q3 verwenden: d uN F3 = −k . d x Q3 4.2.6

Probleml¨osung - Bestimmung der globalen Knotenwerte Q1 , Q2 , Q3

Gleichung (4.31) ist ein lineares Gleichungssystem zur Bestimmung der Unbekannten Q1 , Q2 und F3 (der Knotenwert Q3 ist ja bereits bekannt)! F¨ur die Bestimmung der L¨osung linearer Gleichungsysteme stehen heutzutage verschiedene Verfahren zur Auswahl. F¨ur mehr als 500.000 Unbekannte haben die sog. direkten L¨oser wie das Gauß-Verfahren oder die LU-Dekomposition jedoch zunehmend Schwierigkeiten, so dass man f¨ur große Probleme im Allgemeinen sog. iterative L¨oser, die schrittweise eine L¨osung verbessern, verwendet. Einfache Beispiele daf¨ur sind das Jacobiund Gauß-Seidel-Verfahren, praktisch werden aber eher Mehrgitterverfahren oder vorkonditionierte KrylowUnterraum-Verfahren, wie das Verfahren der konjugierten Gradienten oder GMRES, verwendet17 . Aufgrund der Gr¨oße der Gleichungssysteme ist manchmal der Einsatz von Parallelrechnern n¨otig. In unserem Fall kann man die L¨osung noch per Hand berechnen und wir erhalten

4.2.7

fn L s L2 + + Tb 2k k 3 s L2 fn L + + Tb 8k 2k

Q1

=

Q2

=

F3

=

s L + fn

[ Q3

=

Tb ]

Finite-Elemente-L¨osung

Um die FE-L¨osung auf dem betrachteten Gebiet zu erhalten, m¨ussen die gerade bestimmten Knotenwerte Q1 , Q2 , Q3 in den urspr¨unglichen FE-Ansatz Gl. (4.22) u (x) = N

M X e=1

{N} · {Q} = T

e

[ e

( {n1 n2 }

Q1 Q2

)!

( , {n1 n2 }

Q2 Q3

)! (4.32)

eingesetzt werden. Beachten Sie: • Konstruktionsbedingt ist die L¨osung Gl. (4.32) eine st¨uckweise zusammengesetzte Funktion, die aber trotzdem stetig bzgl. x ist. S ⇒ Daher die Verwendung von e = Vereinigungsmenge nicht-¨uberlappender Elemente. • Die Knotenwerte {Q} sind global durchnummeriert, die FE-Ansatzfunktionen {N} = {n1 n2 } sind jedoch lokal je Element definiert (siehe Gl. (4.21))! 17 Auf einem 1-GFLOP/s-Rechner ben¨ otigt z.B. ein Mehrgitterverfahren zur L¨osung eines Gleichungssystems mit 1 · 107 Unbekannten eine Sekunde; das Gauß-Verfahren w¨urde hierzu 21.000 Jahre ben¨otigen!

Angewandte Mathematik II

39

– Man kann keine geschlossene Formel f¨ur die L¨osung im gesamten Gebiet Ω angeben! – Will man die L¨osung graphisch darstellen, so muss man Gl. (4.32) elementweise zeichnen F¨ur unser Beispiel k¨onnen wir nun den elementweisen Verlauf bestimmen (siehe Abb. 18): Element e = 1: ( {n1 n2 }

Q1 Q2

) =

Element e = 2:

n j (x) Qej

j=1

= mit: XL = 0 ; XR =

2 X

L : ⇔ 2

1 1 (XR − x) Q1 + (x − XL ) Q2 ; XL ≤ x ≤ XR le le  2 L 2 L − x Q1 + xQ2 ; 0≤x≤ L 2 L 2

2 2 L (L − x) Q2 + x− Q3 ; L L 2

L ≤x≤L 2

N

u (x)

Gerade < > lineare Interpolation

0 Q1 1 0 1

1 0 0 1 0 1

Q2 11 00 00 11 00 11

Q3

L 2

0

xL

Ω e=1

xR xL

L Ωe=2

x

xR

Abbildung 18: Beispiel: Lineare Interpolation mittels zweier Geradenst¨ucke Hierzu noch einige Bemerkungen: • Die exakte analytische L¨osung ist quadratisch in x! (siehe Gl. (4.8)) – Die FE-N¨aherung ist st¨uckweise linear! • An den Knoten stimmt die FE-L¨osung mit der analytischen L¨osung u¨ berein! Dazwischen gibt es einen Fehler aufgrund der (vorher festgelegten) linearen Interpolation! • Der W¨armefluss Element-Mitte!)

du dx

an den Knoten stimmt nicht mit der analytischen L¨osung u¨ berein (nur in der

ABER: Am rechten und linken Rand werden die Fl¨usse sogar exakt bestimmt! → GalerkinFormulierung ist dort exakt!

Angewandte Mathematik II

40

• Anhand von Abb. 18 erkennt man sofort, dass die FE-N¨aherung insgesamt verbessert werden kann, wenn man mehr Elemente verwendet. Im Hinblick auf die oben durchgef¨uhrten Schritte zu Entwicklung der FH-L¨osung ist das aber ganz einfach! Erweiterung des FE-L¨osung von M = 2 auf beliebig viele Elemente: 1. Ausgangspunkt ist Gl. (4.27), die bereits allgemein f¨ur M Elemente (und damit M + 1 Knoten) formuliert worden ist! 2. Setze le =

L M

als elementspezifische Breite f¨ur jedes der M Elemente.

3. Aufstellen des globalen Knoten-Vektors {Q} = (Q1 , Q2 , . . . , Q M+1 )T . 4. Erweitern der 2 × 2-Matrizen zu (M + 1) × (M + 1)-Matrizen durch Auff¨ullen der nicht am Element beteiligten Knoten-Spalten und -Zeilen mit Nullen. 5. Assemblierung der Elementgleichungen zum globalen Gleichungssystem f¨ur den Vektor {Q}. 6. Einsetzen der Randbedingungen. 7. L¨osung des Gleichungssystems f¨ur {Q} (Seminaraufgabe!). Dies ist ein beachtliches Ergebnis und zeigt die Systematik hinter der FEM. Der oben stehende Algorithmus ist allgemein und kann in Form einer Schleife sofort f¨ur beliebig viele Elemente programmiert werden. Auf diese Weise l¨aßt sich die N¨aherungsl¨osung relativ leicht verbessern. Alternativ hierzu gibt es aber in der oben dargestellten Herleitung der FH-L¨osung noch eine weitere M¨oglichkeit Einfluß auf die Qualit¨at der N¨aherung zu nehmen. Wir hatten uns gleich zu Beginn auf lineare FEAnsatzfunktionen festgelegt. Dies ist zwar als Zugang, um die FEM zu verstehen, sehr gut geeignet. Wir begehen hierbei jedoch sehenden Auges“ einen fundamentalen Fehler, da wir damit jede Art von Kurven” ¨ krummung von vornherein u¨ berhaupt nicht erfassen k¨onnen, was aber notwendig ist, wenn man sich den Vergleich der FE-L¨osung unseres Stabbeispiels mit der exakten analytischen L¨osung anschaut! Dar¨uber hinaus gibt es noch mathematisch und praktisch viel gewichtigere Gr¨unde: Lineare Elemente garantieren zwar die Stetigkeit der L¨osung an den Elementknoten, nicht aber die Stetigkeit der Ableitung! H¨aufig sind die Stetigkeitsanforderungen aber h¨oher, z.B. um auch den Fluss du/dx einer Gr¨oße berechnen zu k¨onnen! Dies erfordert aber, dass auch die erste und zweite Ableitung stetig sind! Elemente mit Ansatzfunktionen, die den Stetigkeitsanforderungen der L¨osung gen¨ugen, nennt man in FEM-Kreisen konform. Die Erweiterung der FEM auf Ansatzfunktionen in Form von Polynomen h¨oherer Ordnung ist problemlos m¨oglich wie der n¨achste Abschnitt zeigt.

4.3

FE Ansatzfunktionen - Formfunktionen (engl. shape functions)

Der Begriff Formfunktion bezieht sich auf die Form (mathematisch + optisch) der Ansatzfunktionen auf den Elementen. Man verwendet beispielsweise den (saloppen) Begriff lineare Elemente“ f¨ur Elemente, auf denen ” lineare Formfunktionen, wie in unserer Herleitung oben, definiert sind. Die mathematische Form der Ansatzfunktionen ergibt sich aus der Bedingung, dass I: die gesuchte N¨aherungsl¨osung uN (x) z.B. linear oder quadratisch zwischen den Knoten des Elements verlaufen soll und II: an den Knoten des Elements die Werte {Q}e annimmt. Wir wollen dies im Folgenden beispielhaft an linearen“ und quadratischen“ Elementen demonstrieren. ” ”

Angewandte Mathematik II

41

4.3.1

Lineare Formfunktionen - Lineare Elemente“ ” Mathematisch formuliert ergeben die o.g. Bedingungen je Element ueN (x) =

I: II:

α1 + α2 x

=

ueN (xL )

und

Q1

und

(4.33)

ueN (xR )

= Q2

und damit Q1

=

α1 + α2 x L

Q2

=

α1 + α2 x R

Aufl¨osen nach α1 und α2 und Einsetzen in I: mit le = xR − XL ergibt:

=

1 1 (xR − x) Q1 + (x − xL )Q2 le le n1 (x) Q1 + n2 (x) Q2

=

{N}T {Q}e

ueN (x) =

Wir erhalten die bereits bekannten FE-Ansatzfunktionen {N} Gl. (4.21) (siehe Abb.19)!

uN(x)

{N} u 2(x)

u 1(x)

Q2

1 11 00 00 11 00 11

xL

11 00 00 11 00 11

uN(x)

Q1

Formfunktion im Element Ωe

11 00 00 11 00 11

11 00 00 11 00 11

xR x

Ωe

xR x

xL N

FE−Darstellung von u (x) im Element Ωe

Abbildung 19: Lineare Ansatzfunktionen und die darauf basierende Darstellung der FE-N¨aherungsl¨osung uN (x) im Element Ωe .

4.3.2

Quadratische Elemente

Die Vorgehensweise zur Bestimmung quadratischer Ansatzfunktionen ist exakt dieselbe wie bei linearen Elementen, mit dem Unterschied dass man jetzt eine allgemeine quadratische Funktion zugrunde legt, also I:

ueN (x) = α1 + α2 x + α3 x2

Das Problem ist allerdings jetzt, dass man drei Parameter α1 , α2 , α3 bestimmen muss, daf¨ur aber zwei Knoten Q1 , Q2 nicht ausreichen! Die L¨osung f¨ur diese Problem ist aber ganz einfach! Man einfach noch einen weiteren Knoten hinzu, z.B. in der Mitte des Elements! Damit ergibt sich I: II:

ueN (x) = ueN (xL ) =

α1 + α2 x + α3 x2 und le Q1 , ueN (xL + ) = Q2 2

(4.34) und ueN (xR ) = Q3

Angewandte Mathematik II

42

und damit Q1

=

Q2

=

Q3

=

α1 + α2 xL + α3 x2L le le α1 + α2 (xL + ) + α3 (xL + )2 2 2 α1 + α2 xR + α3 xR2

Grunds¨atzlich kann diese Gleichungen nach den drei Parametern α1 , α2 , α3 aufl¨osen, man macht sich dabei das Leben allerdings unn¨otig schwer! Warum? Die Wahl der Koordinatenbezeichnungen f¨ur die R¨ander eines Elements (XL und XR ) war zwar anschaulich als wir die FEM hergeleitet haben, ist aber f¨ur praktische Berechnungen umst¨andlich. Beim Betrachten von Gl. (4.23) – quasi dem Kernst¨uck der FEM – f¨allt auf, dass s¨amtliche Berechnungen, die Auswertung der Integrale und der R¨ander, unabh¨angig von einander auf jedem Element durchgef¨uhrt werden. Diese Integrale sind allesamt bestimmte Integrale und die Grenzen sind die R¨ander der Elemente. Man darf sich, wie in solchen F¨allen u¨ blich, den Koordinatenursprung legen wie man will. Dabei w¨ahlt man u¨ blicherweise eine Lage, die die mathematische Formulierung der Ansatzfunktionen m¨oglichst einfach macht! Bevor wir also die α1 , α2 , α3 bestimmen, f¨uhren wir eine simple Koordinatentransformation durch, in dem wir XL ≡ 0 setzen. Damit bekommt der mittlere Knoten die Koordinate le /2 und der rechte Knoten le . Damit vereinfacht sich Gl. (4.35) zu =

α1

Q2

=

le le α1 + α2 + α3 2 2

Q3

=

α1 + α2 le + α3 l22 .

Q1

!2

Aufl¨osen nach α1 , α2 und α3 und Einsetzen in I: ergibt: ! ! ! x2 x x x x x N Q1 + 4 1 − Q2 + 2 − 1 Q3 ue (x) = 1 − 3 + 2 le le le le le le = n1 (x) Q1 + n2 (x) Q2 + n3 (x) Q3 =

{N}T {Q}e

Der Verlauf der drei Ansatzfunktionen ist gemeinsam mit der Koordinatentransformation in Abb.20 dargestellt. Was w¨urde sich an unserem Beispiel a¨ ndern, wenn wir quadratische statt lineare Ansatzfunktionen verwenden? Man kann dies leicht nachvollziehen. Beginnend mit der Berechnung der Ableitungen der Ansatzfunktionen Gl. (4.25) ( ! ! !) d {N} 1 x x x = 4 − 3 ;4 1 − 2 ; 4 −1 (4.35) dx le le le le berechnet man die Integrale (dies ist in dem Fall etwas aufw¨andiger, weil in der Matrix nicht nur Konstanten auftauchen, sondern Funktionen von x) und Randterme. Der einzige Unterschied besteht darin, dass nun je Element 3 × 3 Matrizen auftauchen und f¨ur das gesamte Gebiet die L¨osung f¨ur 5 Knoten (3 Endknoten und 2 Mittelknoten) {Q} berechnet werden muss, so dass bei der Assemblierung die 3 × 3 Matrizen zu 5 × 5 Matrizen aufgebl¨aht“ werden m¨ussen. ” Nach Einsetzen der Randbedingungen kann auch dieses Gleichungssystem gel¨ost werden. Die FE-L¨osung ist wieder eine st¨uckweise Darstellung allerdings jetzt mit je drei quadratischen Ansatzfunktionen pro Element. Im Vergleich mit der analytischen L¨osung stellt man fest, dass beide exakt u¨ bereinstimmen, was allerdings nur Zufall ist und im Allgemeinfall nicht so ist. Fazit: Nach dem selben Schema lassen sich nun auch Ansatzfunktionen 3. und h¨oherer Ordnung definieren. Entsprechend der Anzahl der zu bestimmenden αi m¨ussen entsprechend mehr Knoten (zus¨atzlich zu den Randknoten) im Element hinzugef¨ugt werden, um ausreichend Bedingungen zu bekommen. Schaut man sich die

Angewandte Mathematik II

43

Abbildung 20: Quadratische Ansatzfunktionen im Element Ωe und Koordinatentransformation der lokalen Elementr¨ander.

Konstruktion bzw. die zugeh¨origen Schaubilder Abb.20 und Abb.19 an, so f¨allt auf, dass die so konstruierten Formfunktionen alle die Eigenschaft    1 i = j (4.36) ni (X j ) = δi j und δi j =   0 i , j teilen, wobei X j die Knotenkoordinaten sind. Mit anderen Worten, nur an einem Knoten X j ist der Wert der Ansatzfunktion ni (X j ) eins, an allen anderen ist sie Null! Ein solches Interpolationspolynom nennt man LagrangeInterpolationspolynom. Grunds¨atzlich kann man auch andere Polynome verwenden und es gibt in den verschiedenen FEM-Werkzeugen u¨ blicherweise einen ganzen Zoo“ von Polynomtypen. Je nach Wahl der Ansatz-Funktionen wird die Form“ ” ” der N¨aherungsl¨osung beeinflusst, was entscheidenden Einfluß auf die Qualit¨at der N¨aherung (siehe Abb. 21) sowie auf die Weiterverarbeitung“ (Stichwort: Differenzierbarkeit!) hat. Man muss aber auch beachten, dass mit ”

Abbildung 21: Ab¨angigkeit der FEM-N¨aherungsl¨osung von Typ der Formfunktion.

zunehmender Ordnung auch mehr Knoten je Element berechnet werden m¨ussen und damit der Rechenaufwand steigt18 !

18 Bei

dreidimensionalen Elementen kommen schnell 50 bis 80 Knoten pro Element zusammen!