�

WESTFÄLISCHE WILHELMS-UNIVERSITÄT

MÜNSTER Mathematisches Institut

Numerische Analyse und Simulation des Blutflusses durch künstliche Herzklappen

Sabine Aust

Münster im März 2000

Inhaltsverzeichnis

Einleitung 1

1 Das menschliche Herz 4 1.1 Aufbau und Funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Diastole und Systole . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Herzklappen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Künstliche Herzklappen 12 2.1 Angeborene und erworbene Klappenfehler . . . . . . . . . . . . . . . 12 2.2 Die wichtigsten Prothesetypen . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Schwachstellen der unterschiedlichen Prothesen . . . . . . . . . . . . 17

3 Physikalische Grundlagen 21 3.1 Der Idealfall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Historische Ansätze . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Ein moderner Ansatz: Mary J. King . . . . . . . . . . . . . . . . . . . 24

4 Navier-Stokes-Gleichungen und Reynoldszahl 29

4.1 Herleitung der Euler-Gleichungen . . . . . . . . . . . . . . . . . . . . 29

4.1.1 Beispiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Herleitung der Navier-Stokes-Gleichungen . . . . . . . . . . . . . . . 35 4.3 Die Reynolds-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5 Theoretische Vorbereitungen 41

5.1 Funktionenräume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.2 Variationsgleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5.2.1 Beispiel 1: Lineares Randwertproblem . . . . . . . . . . . . . . 45 5.2.2 Beispiel 2: Poissonsche Differentialgleichung . . . . . . . . . . 46 5.2.3 Ausblick auf Lösungsmethoden . . . . . . . . . . . . . . . . . 47

5.3 Das Variationsproblem im Hilbert-Raum . . . . . . . . . . . . . 48 5.3.1 Satz von Lax-Milgram . . . . . . . . . . . . . . . . . . . . . . 49 5.3.2 Anwendung auf ein Modellproblem . . . . . . . . . . . . . . . 51 5.3.3 Das Optimierungsproblem . . . . . . . . . . . . . . . . . . . . 53

6 Diskretisierung der Probleme 55 6.1 Galerkin-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.2 Erweiterungen auf nichtlineare Probleme . . . . . . . . . . . . . . . . 58 6.3 Methode der finiten Elemente . . . . . . . . . . . . . . . . . . . . . . 62

6.3.1 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.3.2 Methode der gewichteten Residuen . . . . . . . . . . . . . . . 64 6.3.3 Modellproblem: Die 1D Poisson-Gleichung . . . . . . . . . . . 65 6.3.4 Elementmatrizen . . . . . . . . . . . . . . . . . . . . . . . . . 68

6.4 Iterative Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.4.1 Klassische Methoden: Jacobi und Gauss-Seidel Verfahren … 71 6.4.2 Relaxationsverfahren: SOR und SSOR . . . . . . . . . . . . . 74 7 Existenz- und Eindeutigkeitssätze zu den Navier-Stokes Gleichun-

gen 76 7.1 Klassische Resultate: R. Temam, 1977 . . . . . . . . . . . . . . . . . 76

7.1.1 Notationen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.1.2 Schwache Lösungen des Linearen Problems . . . . . . . . . . . 77

7.2 Aktuelle Resultate: M.Wiegner, 1999 . . . . . . . . . . . . . . . . . . 82 7.2.1 Der Stokes-Operator und die Stokes-Halbgruppe . . . . . . . . 82 7.2.2 Schwache Lösungen des Nichtlinearen Problems . . . . . . . . 85 7.2.3 Existenz und Eigenschaften starker Lösungen . . . . . . . . . 88

8 Implementierung unter “Diffpack 3.0” 97 8.1 Geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 8.2 Penalty-Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 8.3 Lösung des nichtlinearen Systemes . . . . . . . . . . . . . . . . . . . . 100

8.3.1 Successive Substitutions . . . . . . . . . . . . . . . . . . . . . 101 8.3.2 Newton-Raphson-Methode . . . . . . . . . . . . . . . . . . . . 101

8.4 Programmaufbau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 8.5 Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Zusammenfassung 112

A Medizinische Fachbegriffe 112

B Programm-Code und Input-Files 114 B.1 NsPenalty1.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 B.2 NsPenalty.cpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 B.3 main.cpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 B.4 Generierung der Super-Elemente-Gitter . . . . . . . . . . . . . . . . . 128 B.5 Input-File für Super-Elemente-Gitter . . . . . . . . . . . . . . . . . . 130 B.6 Input-File für BOX WITH BELL-Geometrie . . . . . . . . . . . . . . . . 131 Literaturverzeichnis 133

Abbildungsverzeichnis

1 Längsschnitt durch das Herz . . . . . . . . . . . . . . . . . . . . . . . 5 2 Phasen der Herztätigkeit . . . . . . . . . . . . . . . . . . . . . . . . . 9 3 Querschnitt durch das Herz . . . . . . . . . . . . . . . . . . . . . . . 10 4 Doppelklappenersatz . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 5 Klappentypen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6 Strömungsprofile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 7 Komplikationen bei SJM-Patienten . . . . . . . . . . . . . . . . . . . 18 8 Komplikationen bei Carbomedics-Patienten . . . . . . . . . . . . . . . 18 9 Überlebensraten der Carbomedics-Patienten . . . . . . . . . . . . . . 19 10 Überlebenraten der SJM-Patienten . . . . . . . . . . . . . . . . . . . 19 11 Strömung um einen 2-D Zylinder . . . . . . . . . . . . . . . . . . . . 25 12 schematische Darstellung einer Doppelflügelklappe . . . . . . . . . . . 27 13 Formfunktionen N_{i} und deren Ableitungen N′_{i} . . . . . . . . . . . . . 67 14 Zusammensetzung der globalen Matrix aus den Elementmatrizen … 69 15 abstrahierte 2-D Herzklappe . . . . . . . . . . . . . . . . . . . . . . . 98 16 FE-Gitter unter Verwendung der Super-Elemente-Methode . . . . . . 99 17 Druck p auf BOX WTH BELL-Geometrie . . . . . . . . . . . . . . . . . . 104 18 Geschwindigkeit v auf BOX WTH BELL-Geometrie . . . . . . . . . . . . 105 19 Druck p auf Superelemente-Geometrie . . . . . . . . . . . . . . . . . . 105 20 Geschwindigkeit v auf Superelemente-Geometrie . . . . . . . . . . . . 106 21 Druck p auf Superelemente-Geometrie . . . . . . . . . . . . . . . . . . 106 22 Geschwindigkeit v auf Superelemente-Geometrie . . . . . . . . . . . . 107 23 Druck p auf Superelemente-Geometrie . . . . . . . . . . . . . . . . . . 107 24 Geschwindigkeit v auf Superelemente-Geometrie . . . . . . . . . . . . 108 25 Druck p auf Superelemente-Geometrie . . . . . . . . . . . . . . . . . . 108 26 Geschwindigkeit v auf Superelemente-Geometrie . . . . . . . . . . . . 109

Einleitung

In den letzten Jahrzehnten haben sich sowohl Hardware als auch numerische Algorithmen so stark weiterentwickelt, daß es Wissenschaftlern und Ingenieuren verschiedenster Fachrichtungen möglich ist, komplexe Vorgänge am Computer zu simulieren. Selbst in den klassischen Naturwissenschaften wie Biologie und Medizin sind numerische Simulationen zu einem essentiellen Werkzeug der Forschung und Weiterentwicklung geworden.

In dieser Arbeit soll nun speziell das medizintechnische Problem der künstlichen Herzklappen betrachtet werden. Zur genauen Analyse und Simulation des Blutflusses durch eine künstliche Herzklappe bedarf es einiger biologischer, physikalischer und mathematischer Erläuterungen:

Das 1. Kapitel befaßt sich ausschließlich mit den biologischen Gegebenheiten des menschlichen Herzens, speziell der Aortenklappe, während das 2. Kapitel dann die unterschiedlichen Formen möglicher Klappenprothesen beschreibt und deren Vorund Nachteile vergleicht. Die medizinischen bzw. biologischen Fachbegriffe, die hierzu nötig sind, werden im Text durch kursive Schreibweise gekennzeichnet und im Anhang genauer erläutert.

Kapitel 3 und 4 schaffen im Anschluß daran die physikalischen Grundlagen, die zum Verständnis der Klappenbewegungen notwendig sind. Es werden die Ursachen der wirbeligen oder sogar turbulenten Blutströmung skizziert und die beschreibenden mathematischen Gleichungen, die Navier-Stokes Gleichungen hergeleitet.

Die komplexe mathematische Theorie der Variationsgleichungen und Funktionenräume, die allgemein hinter partiellen Differentialgleichungen dieser Gestalt steht, wird in Kapitel 5 vorbereitet, so daß in Kapitel 6 auf Lösungsansätze in geeigneten Funktionenräumen eingegangen werden kann. Im Hinblick auf die spätere Implementierung wird auf das Galerkin-Verfahren, die Methode der finiten Element und auf iterative Methoden zur Lösung der resultierenden Gleichungssysteme besonders eingegangen.

Das Ziel des theoretischen Teils dieser Arbeit wird in Kapitel 7 erreicht. Hier werden zunächst zwei Beweise der Existenz schwacher Lösungen des linearen und nichtlinearen Navier-Stokes-Problemes geliefert, um anschließend einen Beweis der Existenz starker Lösungen der vollständigen Navier-Stokes Gleichungen in der Dimension n ≥ 3 zu erläutern. Diesers letzte Existenzbeweis liefert schon anhand von genügend kleinen Anfangswerten eine globale, eindeutige Lösung. Die Berechtigung der gängigen numerischen Verfahren ist hiermit gesichert.

Die Umsetzung dieser biologischen, physikalischen und mathematischen Grundlagen und Erkenntnisse erfolgt in Kapitel 8, welches die Implementierung eines vereinfachten 2-dimensionalen, stationären Modelles zur Simulation des Blutflusses durch eine künstliche Aortenklappe beinhaltet. Unter Verwendung der C++ Bibliothek “Diffpack 3.0” wird der Aufbau eines klassischen Finite-Elemente-Lösers für die Navier-Stokes Gleichungen beschrieben und auf das stationäre Problem mit unterschiedlichen Positionen der Klappenflügel angewendet. Besondere Aufmerksamkeit soll dabei den Auswirkungen verschiedener Öffnungswinkel gewidmet werden. Zum Vergleich der Ergebnisse dienen sowohl die simulierten, als auch die experimentellen Ergebnisse von Mary J. King aus ihrer Doktorarbeit “Computational and experimental studies of flow through a bileaflet mechanical heartvalve” aus dem Jahre 1994 (siehe [23]).

Der Vergleich zeigt deutlich, daß ein 2-dimensionales zeitunabhängiges Modell nicht exakt die Bildung von Turbulenzen und Wirbeln lokalisieren kann, aber dennoch durchaus in der Lage ist, Tendenzen der Blutströmung in der unmittelbaren Umgebung der mechanischen Herzklappe vorherzusagen.

Danksagung

Die ersten, denen ich danken möchte, sind Prof. Colin Cryer und Prof. Michael Wiegner. Prof. Cryer war derjenige, der meine Begeisterung für das interdisziplinäre Arbeiten geweckt und durch das vergangene Jahr begleitet und gefördert hat. Prof. Wiegner brachte mir durch geduldige und wiederholte Hilfestellungen seinen Beweis zu Theorem 7.1 näher, so daß Kapitel 7 in greifbare Nähe rückte. An dieser Stelle sei für das mir und meiner Arbeit entgegengebrachte Interesse gedankt.

Doch ohne die technische Hilfe vieler gutmütiger Menschen wäre diese Diplomarbeit nicht so bald fertiggestellt worden. Allen voran sind da Dr. Frank Wübbeling und Axel Feldmann zu nennen, die mich immer wieder milde lächelnd mit einem Knopfdruck von mir unlösbar erscheinenden Computer-Problemen befreiten, dicht gefaolg von Jutta Lücking, die sich gemeinsam mit mir durch Diffpack kämpfte und über jedes bunte Bild freute. Insgesamt war die Atmosphäre in dem Institut für numerische Mathematik durch den allgemein humorvollen Umgangston mehr als angenehm und ich bin froh, in diese Etage gelandet zu sein.

Begleitet, unterstützt und motiviert haben mich in den letzten Monaten viele gute Freunde, unter denen Michael G. W. Jacob ein besonderer Dank gebührt. Seine allumfassenden Ratschläge und sein tapferes Korrekturlesen, das von Axel Naumann noch unterstützt wurde, waren mir eine große Hilfe. Dann waren da noch viele wichtige Menschen in meinem Umfeld wie Kerstin Ronneberger, Helke Kläning, Sandra Denningmann und die unverwüstliche und doch so charmante Mensa-Runde samt allen gelegentlichen Gästen, um nur einige von denen zu nennen, die ich vermissen werde.

Zu guter Letzt noch ein Danke, das mit Worten kaum noch auszudrücken ist, an Markus Kirchler, der in den vergangenen Monaten liebevoll über manche Stresserscheinung hinweggesehen hat, und an meine Familie, ohne die ich jetzt nicht da wäre, wo ich bin. Ohne Euch hätte ich selbst nicht geglaubt, daß es zu schaffen ist.

1 Das menschliche Herz

Im folgenden sollen Aufbau, Funktion und krankhafte Veränderungen des menschlichen Herzens in ihren Grundzügen beschrieben werden; vergleiche dazu [12], [33] und [29].

1.1 Aufbau und Funktion

Das menschliche Herz ist ein pulsierender Hohlmuskel, der die Bewegung der Blutflüssigkeit in unserem geschlossenen Blutsystem bewirkt. Form und Größe des Herzens entsprechen in etwa der Faustgröße der jeweiligen Person. Seine Oberseite, die Herzbasis, an welcher die großen Gefäße münden bzw. entspringen, hat einen deulich größeren Umfang als die nach unten links gerichtete Herzspitze. Die Herzbasis liegt unmittelbar hinter dem Brustbein, während die Herzspitze nach links von der Mittellinie abweicht.

Der größte Teil des Herzens ist vom Perikard, dem Herzbeutel, eingehüllt. Das Perikard bildet eine allseits geschlossene seröse Höhle mit normalerweise nicht mehr als 20 ml Flüssigkeit im Binnenraum. Wie andere seröse Höhlen auch, dient der Herzbeutel dazu, die Reibung zwischen dem beweglichen Herzen und den mehr oder weniger fixierten Nachbarorganen herabzusetzen (ähnlich zweier Glasplatten mit einem Tropfen Wasser dazwischen). Außerdem überträgt dieser Herzbeutel den Sog des Lungengewebes auf das Herz und sorgt somit für einen Druckausgleich im gesamten Thoraxraum.

Das Herz wird von einer muskulösen Scheidewand, dem Septum, in eine linke und eine rechte Herzhälfte getrennt, welche sich jeweils in Vorhof (Atrium) und Kammer (Ventrikel ) unterteilen. Das Septum ist im Vorhofbereich dünn und sehnig, zwischen den Kammern aber außerordentlich muskelstark, da die mechanische Beanspruchung hier deutlich höher ist.

Der rechte Vorhof erhält durch die Hohlvenen sauerstoffarmes Blut aus dem Körperkreislauf und befördert es in die rechte Kammerhälfte. Von hier aus wird es zum Gasaustausch über die Lungenarterie in die Lunge gepumpt, gelangt nach deren Passage als sauerstoffreiches Blut in die Lungenvenen und dann in den linken Vorhof (”kleiner“ Blutkreislauf). Aus dem linken Vorhof wird es in die linke Kammer, von dieser mit hohem Druck in die Aorta und weiter in den Körperkreislauf gepumpt (”großer“ Blutkreislauf).

Wie das gesamte Gefäßsystem, so hat auch das Herz einen dreischichtigen Aufbau:

Abbildung 1: Längsschnitt durch das Herz (aus [33], S.127) 5 • Das Endokard, die glatte Innenschicht des Herzens, kleidet den Blutraum

vollständig aus und ist somit entscheidend für eine gute Gleitfähigkeit des Blutes. Schon kleine Unebenheiten, z.B. aufgrund einer Endokarditis, erhöhen den Strömungswiderstand relevant. Das Blut fließt an diesen Stellen dann langsamer, es können sich Thrombozyten ablagern und ineinander verkeilen - es enstehen Blutgerinsel, bzw. Thromben. Wenn sich ein Thrombus ablöst und in den allgemeinen Kreislauf gelangt, kann er eine wichtige Arterie, z.B. im Gehirn, verstopfen und damit einen Schlaganfall verursachen (Embolie).

Unterhalb dieses einschichtigen Gewebes finden sich reichlich elastische Fasern. Aus einer besonders kollagenfaserreichen Falte des Endokards bestehen auch die Herzklappen und Sehnenfäden.

• Das Myokard, die Herzmuskelschicht, schließt sich nach außen an das Endokard an. Da das Myokard die eigentliche Herzarbeit verrichten muß, ist es die dickste Schicht der Herzwandung. Es besteht aus einkernigen Herzmuskelzellen, die mechanisch und elektrisch miteinander verbunden sind. Im Ventrikelbereich schrauben sich die Muskelfasern in den äußeren Myokardschichten spiraliförmig von der Herzklappenebene zur Herzspitze hinab, biegen dort um und winden sich innen wieder in Richtung Herzbasis zurück. Dadurch stehen die Muskelfasern der äußeren und inneren Myokardschicht häufig senkrecht aufeinander. Die Aktin- und Myosinfilamente der Herzmuskulatur überlappen sich also in Ruhe so stark, daß sie sich tendenziell gegenseitig behindern. Bei körperlicher Belastung erhöht sich das Volumen der Ventrikel, die Muskulatur wird besser vorgedehnt und arbeitet effektiver, wodurch sich Kraft und Schlagvolumen des Herzens erhöhen.

Diese autonome Anpassung der Muskulatur besitzt eine besonder große Bedeutung bei der Herzinsuffizienz: Ein geschwächter Herzmuskel neigt stets dazu, ein größeres Volumen anzunehmen (”ausgelatschtes Herz“), was den oben beschriebenen Mechanismus aktiviert und die Herzschwäche über lange Zeit hinweg ausgleichen kann.

• Das Epikard sorgt für eine gute Gleitfähigkeit und Beweglichkeit des Herzens in seiner Umgebung. Es stellt gleichzeitig die äußerste Schicht des Herzens und das innere Blatt des Herzbeutels dar.

Wegen der extremen Dicke der Herzwand kann das Herz nicht direkt aus dem durchgepumpten Blut versorgt werden — die Diffusionsstrecken wären viel zu lang. Es ist ein eigenes Blutversorgungssystem nötig, welches von den Herzkranzgefäße gebildet wird. Bei Verschluß einer Herzkranzarterie oder eines größeren Arterienastes wird die von diesem Gefäß versorgte Muskulatur nicht mehr mit Sauerstoff und Nährstoffen beschickt und stirbt ab. Auf diese Weise ensteht ein Herzinfarkt.

1.2 Diastole und Systole

Um eine Flüssigkeit durch eine Leitung zu pumpen, gibt es zwei sich grundlegend unterscheidende Möglichkeiten: die Kreiselpumpe und die Ventilpumpe. Da die biologische Evolution das Rad und damit auch drehbare Systeme aller Art niemals entwickelt hat, muß das Herz eine Ventilpumpe sein, die ähnlich wie ein Blasebalg funktioniert: Ein- und Auslaßventil entsprechen den Herzklappen, der eigentliche Blasebalg ist die Herzkammer.

Bei der rhythmischen Tätigkeit des Herzens unterscheidet man in einem Zyklus zwischen der Kontraktions- und Austreibungssphase (Systole) sowie der Erschlaffungs- und Füllungsphase (Diastole).

1. Systole

• Kontraktionsphase: Die Systole beginnt mit einer Anspannungsphase von ca. 60 ms, in der alle Herzklappen geschlossen sind. Zunächst spannt sich die Ventrikelmuskulatur an und übt einen Druck auf das Blutvolumen in den Ventrikeln aus. Das Volumen der Ventrikel ändert sich dabei nicht, da das Blut wie jede wässrige Flüssigkeit praktisch inkompressibel ist. Daher wird die Anspannungsphase auch isovolumetrische Kontraktion genannt. Der Ventrikeldruck nimmt so lange zu, bis er den Druck in der Aorta bzw. in dem Truncus pulmonalis übersteigt. Der rechte Ventrikel hat es dabei deutlich leichter, da er nur den diastolischen Pulmonalarteriendruck von 10 mmHg (Aorta 80 mmHg) überwinden muß. Die Anspannungsphase endet daher rechts früher als links. Zu Beginn der Anspannungsphase machen beide Ventrikel eine Formveränderung durch: Aus der länglichen Gestalt während der Diastole wird eine annähernde Kugelform. Diese Formveränderung teilt sich der Brustwand als dumpfer Schlag mit und ist als erster Herzton zu hören und über der Herzspitze zu tasten.

• Austreibungsphase: Nachdem der Ventrikulardruck die diastolischen Werte in Truncus pulmonalis und Aorta überschritten hat, öffnen sich die Klappen zu den beiden großen Arterien und es beginnt die ca. 200 ms lange Austreibungsphase, also die eigentliche Bewegung der Blutsäule in den großen Arterien. Die Einlaßventile zwischen Ventrikel und Atrium bleiben wegen des hohen Drucks im Ventrikel weiterhin geschlossen. Trotz des ständigen Auswurfes von Blutvolumen steigt der Ventrikeldruck während der ersten zwei Drittel der Austreibungsphase weiterhin an und erreicht Spitzenwerte von ca. 130 mmHg (links) bzw. 25 mmHg (rechts).

Normalerweise treiben beide Ventrikel nur knapp die Hälfte ihres Inhaltes in die großen Arterien aus; das Schlagvolumen beträgt durchschnittlich 70 ml, kann bei körperlicher Belastung jedoch deutlich gesteigert werden. Durch den Blutauswurf verkleinert sich natürlich das Ventrikelvolumen. Dieses verringerte Volumen zieht einerseits die Herzspitze geringfügig in Richtung Herzbasis, andererseits senkt sich die Klappenebene in Richtung Herzspitze. Die Bewegung dehnt und öffnet die Vorhöfe; venöses Blut wird in das rechte und linke Atrium ”eingesogen“ (Ventilmechanismus).

2. Diastole

• Erschlaffungsphase: Die Diastole beginnt mit einer Erschlaffungsphase von ca. 40 ms. Die Muskulatur der Ventrikelwand verliert ihre Spannung, der Ventrikeldruck sinkt unter den Blutdruck in der Aorta bzw. den Pulmonalarterien ab. Nun schließen sich auch die Ausflußklappen wieder; man spricht von isovolumetrischer Entspannung, wobei sich das eigentliche Zuschlagen durch einen hell klingenden zweiten Herzton verrät. Die Ventilebene wandert geringfügig in Richtung Herzbasis, wird aber durch das große Atriumvolumen an einer vollständigen Hebung gehindert.

• Füllungsphase: Ist die Herzmuskulatur vollständig erschlafft, fällt der Druck im Ventrikel unter den im Atrium. Die Einlaßklappen öffnen sich und es beginnt eine Füllungsphase von sehr variabler Dauer, die mit einer erneuten Systole endet. Da die Atrien aufgrund des Ventilebenenmechanismus prall mit Blut gefüllt sind, findet am Beginn der Füllungsphase ein schneller Einstrom in die Ventrikel statt, der mit verschwindendem Druckunterschied allmählich nachläßt. Am Ende der Füllungsphase kontrahiert sich das Atrium und trägt so aktiv zur Vergrößerung des Blutvolumens in der Herzkammer bei.

Die rhythmische Fortbewegung der Blutflüssigkeit bewirkt in den nachfolgenden zentralen Gefäßabschnitten pulsierende Schwankungen des Wand- bzw. Blutdruckes. Deshalb werden allgemein für den Blutdruck zwei Werte, der systolische und der diastolische Druck angegeben. Um die starken Blutdruckschwankungen zwischen Diastole und Systole zu mildern, sind die großen Arterien mit sehr elastischen Wänden ausgestattet, durch die sie als sogenannte Windkessel fungieren:

Der Begriff ”Windkessel“ stammt aus der Dampfmaschinentechnik und bezeichnet eine technische Lösung des Problems, den stoßweisen Ausstrom von Flüssigkeit aus einer Pumpe in eine gleichmäßige Strömung umzuwandeln. Der physiologische Windkessel besteht im Ausdehnen bzw. Zusammenziehen der Wandung von Aorta und großen Arterien. Das Herz wirft in der Systole Blut aus, während der Blutstrom in der Diastole zum Erliegen kommt. Die Windkesselarterien dehnen sich während der Systole aus und nehmen so eine Teil des Schlagvolumens zusätzlich auf, der in der Diastole passiv wieder ausgepresst wird. Je elastischer die Arterien sind und je größer ihr Binnenvolumen im Verhältnis zum Schlagvolumen ist, desto gleichmäßiger strömt das Blut in die Peripherie ab.

Abbildung 2: Phasen der Herztätigkeit (aus [33], S. 131)

In der oberen Hälfte von Abb. 2 erkennt man die Bewegung der Ventilebene in der Austreibungs- und Füllungsphase. Die untere Hälfte zeigt die Druckverhältnisse in der linken Herzkammer (rot) und in der Aorta (grün) während der Herzaktion.

1.3 Herzklappen

Wie zuvor beschrieben, verhält sich das Herz wie eine Ventilpumpe, wobei die Funktion der Ein- und Auslaßventile von der Herzklappen übernommen wird. Je nach Abbildung 3: Querschnitt durch das Herz (aus [33], S. 128)

10 Aufbau und Funktionalität unterscheidet man grundsätzlich zwischen Segelklappen und Taschenklappen.

Zwischen Vorhöfen und Kammern liegen die beiden großen Segelklappen (siehe Abb. 3): Die Mitral- und die Trikuspidalklappe.

• Die zweizipfelige Mitralklappe liegt zwischen dem linken Atrium und dem linken Ventrikel. Sie besteht aus zwei dünnen, von glattem Endothel überzogenen Einzelsegeln: einem größeren Aortensegel an der Vorderwand und einem kleineren hinteren Wandsegel.

• Ihr Gegenstück im rechten Teil des Herzens ist die dreizipfelige Trikuspidalklappe mit je einem vorderen, hinteren und Septumsegel.

Die Segelklappen sind Einwegventile, die einen Blutstrom nur vom Atrium in den Ventrikel, nicht aber in umgekehrter Richtung zulassen. Zwei weitere Klappen, sogenannte Taschenklappen(sieh Abb. 3), trennen die linke bzw. rechte Herzkammer von den großen Arterien (Aorta und Truncus pulmonalis). Sie sorgen dafür, daß zwischen den einzelnen Herzschlägen kein Blut aus den Arterien in die Ventrikel zurückfließen kann.

• Je nach der verschlossenen Abgangsarterie heißt die rechte Taschenklappe Pulmonalklappe,

• die linke dementsprechend Aortenklappe.

Die Segelklappen sind großflächiger als die Taschenklappen und benötigen daher Sehnenfäden, die von ihrem freien Rand in die Herzkammer ragen und ein ”Durchschlagen“ der Klappen während der Systole in die Vorhöfe verhindern. Sie sind an der Ventrikelinnenwand über Papillarmuskeln befestigt, welche die Länge der Sehnenfäden dem Kontraktionszustand der Kammern anpassen können. Die Taschenklappen bestehen jeweils aus drei identischen Teilen, die sich wie Schöpfkellen in Richtung der Herzkammern vorwölben. Die Klappen sind zwischen zwei Herzschlägen, in der Diastole, durch den Druckunterschied zwischen Arterien und Herzkammern geschlossen; sie stoßen in der Form eines ”Mercedessterns“ aufeinander. Spannt sich der Herzmuskel an (Systole), so drückt das aus dem Herz ausströmende Blut die Taschenklappen an die Wand von Aorta bzw. Truncus pulmonalis; es entsteht eine dreieckige Öffnung. Da die Taschenklappen nur eine relativ kleine Fläche zu verschließen haben, benötigen sie keine Sehnenfäden. Im weiteren wird nur noch die Aortenklappe, also die Taschenklappe zwischen Aorta und linkem Ventrikel betrachtet, da diese das wesentliche Ventil zwischen Herz und Körperkreislauf darstellt.

2 Künstliche Herzklappen

Herzklappen können von Geburt an falsch angelegt oder durch Krankheiten (z.B. Infektionen) geschädigt sein (vergl. [33], S. 129). Weil dadurch der reguläre Blutfluß gestört ist, kann das Blut unter Umständen in die falsche Richtung strömen. Um den Kreislauf aufrecht zu erhalten, muß das Herz deshalb stärker arbeiten.

Diese Situation kann zum Herzversagen führen. Die geschädigte Herzklappe muß daher entweder operativ behandelt oder durch eine neue ersetzt werden. Ein chirurgischer Eingriff am Herzen löst bei den meisten Patienten allerdings große Ängste und Unsicherheit aus. Die veränderten Lebensbedingungen nach der Operation können diese Empfindungen noch verstärken.

So wird von den Betroffenen — besonders in der ersten Zeit nach der Operation — das ”Klicken“ der künstlichen Herzklappen sehr bewußt wahrgenommen. Je nach Herzklappe (Kunststoffklappe oder biologische Klappe) kann es zudem erforderlich sein, daß die Patienten in ihrem weiteren Leben regelmäßig Medikamente einnehmen.

2.1 Angeborene und erworbene Klappenfehler

Zu den seltenen, angeborenen Klappenfehlern (vergl. [21], S. 127-156) zählen die sogenannten Stenosen, wobei am häufigsten noch die angeborenen Aortenklappenstenosen vorkommen. Es handelt sich hierbei um eine Entwicklungsstörung der Klappenflügel. Statt mit drei Flügeln ist die Aortenklappe bei diesen Patienten durch einen genetischen Fehler nur bikuspidal ausgebildet.

Neben den angeborenen Stenosen ist als Hauptursache von Herzklappenerkrankungen das rheumatische Fieber zu nennen. Die rheumatische Endokarditis führt zu einer Verklebung der Klappenflügelränder, woraus nicht selten eine bikuspidale Klappe entstehen kann. Daraufhin kommt es in der Regel zu einer Kalkeinlagerung durch Wirbelbildungen und zu einer hohen Blut-Flußgeschwindigkeit an der Klappe. Die Öffnung der Aortenklappe verengt sich noch mehr. Die analoge Krankheitsentwicklung kann man bei angeborenen Aortenstenosen beobachten.

Je weiter die Klappenstenose ausgebildet ist, desto größer wird das Druckgefälle zwischen linkem Ventrikel und Aorta während der Systole. Der erschwerte Ausfluß des Blutes aus dem Ventrikel in die Aorta führt also zu erhöhtem Druck und somit zu einer starken Hypertrophie der Ventrikelmuskulatur, die dadurch einen erhöhten Sauerstoffbedarf hat.

Im weiteren Verlauf der Erkrankung sinken Herzminutenvolumen und Schlagvolumen in Ruhe ab, wodurch der Weg zur Herzinsuffizienz vorgezeichnet ist. Um einen frühzeitigen Tod durch Myocardversagen, Herzrhythmusstörungen oder Herzinfarkt zu vermeiden, läßt sich die Implantation eines Herzklappenersatzes kaum umgehen.

2.2 Die wichtigsten künstlichen Herzklappentypen und ihre charakteristischen Eigenschaften

Sowohl Aortenklappe als auch Mitralklappe können durch eine künstliche Prothese operativ ersetzt werden, was die Abbildung 4 zeigt:

Abbildung 4: Doppelklappenersatz (aus [6])

Allgemein unterscheidet man zwischen mechanischen und biologischen Prothesen (vergl. [25], S. 25ff und [6]). Jede dieser Prothesen zeigt charakteristische Vorund Nachteile. Zu den mechanischen Prothesen gehören Ball-, Hubscheiben-, Kippscheiben- und Doppelflügelprothesen.

• Bei den Ball- und Hubscheibenprothesen (siehe Abb. 5, a und b) bewegt sich eine Kugel bzw. eine Scheibe in der Öffnungsphase in der Richtung des Abbildung 5: Mechanische Prothesen und Bioprothesen (aus [25], S.26)

14 Blutflusses und wird von einem Käfig festgehalten. In der Schlußphase liegt der Schließkörper auf dem Prothesenring auf. Das Flußprofil der Kugel-Käfig-Prothesen weist einen exzentrischen Verlauf auf, da der Blutstrom von der sich verlagernden Kugel seitlich zu den Wänden abgedrängt wird. Das diastolische Flußmuster einer Hubscheibenklappe zeigt ein ähnliches Flußprofil. Der diastolische Einstrom wird an der Scheibe abgelenkt, seitlich der Scheibe fließt das Blut entlang der Wände nach apikal. Die kurze Flußumkehr in Richtig des Herzens zum Schließen der Klappe erfolgt dann überwiegend zentral auf die Scheibe zu. (siehe Abb. 6, b und c)

• Die Kippscheibenprothesen (siehe Abb. 5, c) besitzen eine bewegliche Scheibe auf einem exzentrisch gelegenen Drehschanier. Durch diese Konstruktion können die meisten Kippscheibenklappen nur bis ca. 70^{∘} geöffnet werden, so daß der Ein- und Ausstrom abhängig von der Einbaurichtung der Prothese unterschliedliche Flußmuster aufweist.

Man unterscheidet einen stärkeren Hauptstrom zur einen und einen schwächeren Nebenstrom zur anderen Wand. (siehe Abb. 6, d und e)

• Ein anderes mechanisches Modell ist die Doppelflügelprothese (siehe Abb. 5, d). Hierbei öffnen sich zwei Klappendeckel wie Türflügel mit einem fast rechten Winkel, ca. 85^{∘}, so daß das Strombahnhindernis im inneren Prothesenbereich nicht sehr ausgeprägt ist. Wenn bei maximaler Öffnung die Flügel beinahe senkrecht stehen, kann das Blut fast ungehindert durch zwei breitere seitliche und einen schmalere zentrale Öffnung auströmen. Die Flußumkehr erfolgt entlang der Seitenwände. (siehe Abb. 6, f und g)

• Die implantierten Bioprothesen (siehe Abb. 5, e und f) bestehen meist aus Schweineaortenklappen oder Kalbsperikard und werden in der Regel auf flexible Bügel aufgenäht, wobei die Basis der Klappe von einem Metallring umgeben wird. Das biologische Gewebe ist an einem Nahtring mit Haltestreben befestigt und zeigt im intakten Zustand fast denselben Bewegungsablauf wie eine normale, menschliche Herzklappe. (siehe Abb. 6, h bis m)

Die weiteren Ausführungen werden sich von jetzt an nur noch auf Doppelflügelprothesen beziehen, da diese in der modernen Praxis am häufigsten eingesetzt werden. Obwohl die Bioprothesen haemodynamisch weitaus günstiger sind, haben sie doch im Vergleich zu den Doppelflügelprothesen einen entscheidenden Nachteil: Wie jedes biologisches Gewebe sind sie nicht unbegrenzt haltbar. Nach ca. 10-15 Jahren zeigen Bioprothesen die ersten Verkalkungs- und Ermüdungserscheinungen, so daß sie in der Regel nur bei Patienten über 60 Jahre eingesetzt werden. Bei jüngeren Patienten wäre das Risiko von weiteren Operationen zu hoch. Da es das Ziel ist, das natürliche Abbildung 6: Strömungsprofile verschiedener Klappenprothesen (aus [25], S.28)

16 Flußprofil so gut wie möglich zu immitieren, erscheinen alle Varianten bis auf die Doppelflügelprothesen und Bioprothesen als nicht besonders geeignet. Weicht das Flußprofil zu sehr vom natürlichen Fall ab, so sind die Belastungen für manche Teile der Aortenwand zu hoch und es entstehen Stellen zu hoher und Stellen zu geringer Flußgeschwindigkeit.

Die gängigsten Modelle der Doppelflügelprothesen sind die St. Jude Medical (kurz: SJM) und die Carbomedics Herzklappen, die sich im Wesentlichen nur durch leicht verschiedene Öffnungswinkel unterscheiden:

Die Carbomedics-Klappen zeigen in maximaler Position einen Öffnungswinkel von 78^{∘}, die St. Jude-Klappen sogar 84^{∘}. Wie sich später heraus stellen wird, rufen diese auf den ersten Blick geringfügigen Unterschiede deutliche Abweichungen in den Flußund Geschwindigkeitsprofilen hervor.

2.3 Schwachstellen bzw. krankhafte Veränderungen des Herzens durch Klappenprothesen

Die häufigsten und wohl auch schwerwiegendsten Krankheitsbilder aufgrund von Herzklappenoperationen sind Hirnblutungen, Magen-Darm-Blutungen, Thromboembolien, haemolytische Anaemien und eine evtl. noch verbleibende Herzinsuffizienz.

• Hirn- und Magen-Darm-Blutungen können als Nebenwirkung von starken gerinnungshemmenden Medikamenten enstehen, die aber meist verabreicht werden müssen, um die erhöhte Trombosegefahr in der Umgebung der Klappenprothese zu kontrollieren.

• Das andere Extrem findet man bei Patienten mit Thromboembolien:

In der Umgebung der Klappenprothese bilden sich Blutgerinsel, die auseinanderbrechen und vom Blutstrom weitergetragen werden. Diese ”Bruchstücke“ können sich festsetzen und z.B. eines der lebensnotwendigen Herzgefäße verstopfen.

• Durch verstärkte Haemolyse an den Herzklappenflügeln kann der Anteil der roten Blutkörperchen stark abfallen. Man spricht von einer haemolytischen Anaemie, der Blutarmut, die eine Störung des Sauerstoff-Transportes durch das Blut verursacht. Alle sauerstoffabhängigen Leistungen, z.B. auch alle Muskel-Leistungen werden auf diese Weise erschwert.

• Sowohl durch die verminderte Sauerstoffversorgung, als auch durch die Thromboseneigung ist das Herz ungewöhnlich großen Belastungen ausgesetzt, so daß es nicht mehr in der Lage ist, das notwendige Blutvolumen in die weiterführenden Gefäße zu pumpen, was zu einer Herzinsuffizienz führen kann.

Abbildung 7: SJM-Herzklappenpatienten ohne spätere Komplikationen, in Prozent; AVR bzw. MVR entsprechen aortic bzw. mitral valve replacement (aus [19], S. 10)

Abbildung 8: links: Carbomedics-Herzklappenpatienten ohne spätere Thromboembolien und ohne spätere Blutungen durch gerinnungshemmende Medikamente (aus [1], S. 631)

Wie man in Langzeitstudien an den St. Jude Medical (siehe Abb. 7 und 10) und den Carbomedics Herzklappen (siehe Abb. 8 und 9) feststellen konnte, findet man in allen Testgruppen nach einigen Jahren Patienten mit den oben beschriebenen Abbildung 9: Überlebensraten der Carbomedics-Herzklappenpatienten nach Aortenklappen-Austausch und nach Mitralklappen-Austausch (aus [1], S. 630)

Abbildung 10: Überlebensraten der SJM-Herzklappenpatienten; AVR bzw. MVR entsprechen aortic bzw. mitral valve replacement (aus [19], S. 10)

Kranheitsbildern, die in vielen Fällen zum Tod führen. Obwohl diese Ereignisse zeitlich stark verzögert auftreten, sind sie immer noch auf die Klappenoperation bzw. auf die klassischen Schwierigkeiten mit der künstlichen Prothese zurückführbar.

Die Mortalitäts-Kurven der beiden betrachteten Klappentypen ähneln sich sehr stark, wobei die Kurve bzgl. der Carbomedics-Prothesen die frühen, postoperativen Todesfälle nicht miteinbezieht, was in der SJM-Graphik der Fall ist. In allen betrachteten Testgruppen fallen diese postoperativen, frühen Fälle mit ca. 7% ins Gewicht.

Der Überblick über diese statistischen Gegebenheiten zeigt, wie häufig doch die unterschiedlichen Krankheitsbilder in der heutigen Zeit noch auftreten. Diese ließen sich theoretisch vermeiden, wenn es möglich wäre, ein naturgetreues Modell einer Herzklappe zu konstruieren, welches keine Unregelmäßigkeiten im Flußprofil mehr aufzeigt. Ein zentrales Kriterium dafür sind die exakte Modellierung und die anschließende Simulation des Blutflusses durch eine mechanische Herzklappe.

3 Physikalische Grundlagen

3.1 Der Idealfall

Betrachtet man die in Abschnitt 1 und 2 beschrieben biologischen und medizinischen Gegebenheiten, so zeichnet sich ab, welche Eigenschaften eine ideale künstliche Herzklappe haben sollte, um dem natürlichen Fall möglichst nahe zu kommen und möglichst wenige Komplikationen zu verursachen (vergl. [23], S. 5):

• Reibungsfreies Öffnen und Schließen der Klappe bereits bei minimaler Druckänderung.

• Keine Druckunterschiede an der Klappenoberfläche während des Blutausstromes.

• Kein Blutrückfluß, während die Klappe sich schließt bzw. ganz geschlossen ist.

• Lebenslange Haltbarkeit ohne Errosionserscheinungen.

• Keine durch die Herzklappe verursachten krankhaften Veränderungen des umliegenden Gewebes.

• Keine Verursachung von Haemolysen und Thrombosen.

Die letzten beiden Bedingungen zielen direkt auf das Flußverhalten des Blutes beim Strom durch die Klappe. Wird die Strömung beim Ausfluß zu turbulent, bilden sich hinter der Klappe, wie wir später noch genauer erkennen werden, starke Wirbel. Diese lassen in den Klappentaschen und in der Aorta einerseits Stellen mit zu hoher und andererseits Stellen mit zu geringer bzw. gar keiner Flußgeschwindigkeit entstehen. Bei zu geringer Flußgeschwindigkeit besteht eine extrem hohe Gerinnungsgefahr, da das Blut fast steht. Zu hohe Geschwindigkeiten lösen meist vermehrte Haemolyse (durch Zerquetschen und Reißen von Blutzellen an den starren Klappenflügeln) und unkontrolliertes Wuchern der arteriellen Wand (durch zu starken Aufprall des Blutstromes) aus.

Um eine ideale Klappenprothese zu konstruieren, braucht man also detaillierte Kenntnisse über das Flußverhalten des Blutes, bzw. Druck und Geschwindigkeit innerhalb des Blutflusses.

3.2 Historische Ansätze

Der Bewegung von Herzklappen widmen sich schon seit Jahrzehnten Physiologen, Physiker und nicht zuletzt Mathematiker. Im Jahre 1912 erkannten Yandell Henderson und F. Elmer Johnson (siehe [18]) als erste die fluiddynamischen Vehältnisse, die zum Öffnen und Schließen der Herzklappen führen und machten dies an einigen erstaunlich einfachen in vitro Experimenten deutlich.

• Die erste Versuchsanordnung besteht aus einem Glasröhrchen, das die Verbindung zwischen einem erhöhten Tintenreservoir und einem darunter liegenden Wasserbehälter bildet. Wird der Tintenzulauf von oben abrupt unterbrochen, kann man beobachten, daß die Tintensäule innerhalb des Röhrchens sofort stillsteht, während der sich schon im Wasser befindende Reststrahl seine Vorwärtsbewegung beibehält.

Dieses produziert in einem kleinen Bereich an der unteren Rohröffnung kurzzeitig einen Unterdruck, so daß klares Wasser von den Seiten in diesen Bereich einströmt.

• Im zweiten Versuch wird zusätzlich noch ein D-förmiges Teilstück in die Mitte des Rohres eingefügt und am unteren Übergang mit einer freischwingenden Membranklappe versehen. Läuft nun zuerst der Tintenfluss gerade das Rohr hinunter in den Wasserbehälter, so verursacht er fast keine Bewegung im gebogenen Teilstück.

Wird dann aber der Fluss z.B. durch kurzes Verschließen des Röhrchens wiederum unterbrochen, so beginnt die Tinte in dem D-Stück zu zirkulieren und veranlaßt die Membranklappe, sich entegegen dem Vorwärtsstrom nach oben zu schließen.

• Der dritte Versuch simuliert nun mit Hilfe einer angepaßten, flexiblen Gummimanschette an der unteren Öffnung des Röhrchnes das Verhalten der Klappenflügel. Die Flüssigkeitssäule innerhalb des Rohres wird deutlich über den Wasserstand des Behälters angehoben und mit einer Deckplatte oben gehalten. Entfernt man die Deckplatte und läßt die Säule nach unten fallen bzw. fließen, so kann man folgendes beobachten:

Fällt der Level innerhalb des Röhrchens unter den des Wasserbehälters, wird der Druck von außen zu groß, so daß die Manschette dann zusammen fällt und das Röhrchen so verschließt.

Mit Hilfe von angefärbter Flüssigkeit und photographischen Momentaufnahmen zeigten Henderson und Johnson recht deutlich, daß die Klappenflügel sich in dem Moment, in dem der Blutstrom unterbrochen wird, nicht wie ”schwingende Türen“ (vergl. [18], S. 81) schließen, sondern sich durch deutliche Wirbel in den Klappenhöhlen von der Basis her ”einrollen“ (vergl. [18], S. 81), während sich das Blut noch im Vorwärtsfluß befindet. Durch diesen Mechanismus wird ein Blutrückstrom in die linke Herzkammer also fast vollständig vermieden. Leider wurden die Schlußfolgerungen von Henderson und Johnson viele Jahre nicht beachtet oder fehlinterpretiert.

Erst im Jahre 1969 veröffentlichten B. J. Bellhouse und F. H. Bellhouse (siehe [2] und [3]) neue fluiddynamische Erkentnisse zur Klappenbewegung, deren Blickwinkel sich etwas von dem von Henderson und Johnson unterschied. Statt, wie diese, die Verantwortung für den Schluß der Klappen ausschließlich dem Effekt eines abbrechenden Flüssigkeitsstrahles zuzuschreiben (”Breaking-of-a-jet-Theorie“), sahen sie die Ursache in der Bildung von Wirbeln an den Klappenflügeln.

• Der Versuchsaufbau sicherte einen laminaren, pulsierenden Fluß durch eine Modellklappe, in dem Wirbelbildung, Geschwindigkeit und Klappenöffnungsgrad zu unterschiedlichen Stadien gemessen und als zeitliche Funktionen aufgezeichnet wurden.

• Die gemessenen Werte lieferten folgende Ergebnisse: Die Stromlinien folgen dem Verlauf der Klappenflügel bis zum distalen Ende des Sinus und teilen sich dann in Sinus-Wirbel (rückwärts) und Hauptstrom (vorwärts). Die Wirbelbildung beginnt in früher Systole, erreicht ihr Maximum am Scheitelpunkt des Herzzyklus und hält bis weit in die Diastole an. Das deutliche Druckgefälle während der ausklingenden Diastole führt zu einer Verlangsamung der Flußgeschwindigkeit in der Aorta, die eine Ausweitung der Sinuswirbel erlaubt, so daß die Klappe sich schließt.

• Erstmals wurden die gemessenen Werte nun auch mit berechneten, also einem mathematischen Modell, verglichen. Auf der Grundlage der physikalischen Gegebenheiten wurden Gleichungen für den Druck aufgestellt, die folgende Werte miteinbezogen: Druckunterschied p_{1} − p_{2}, Geschwindigkeit u, Beschleunigung \frac{∂u}{∂t} (d. h. zeitliche Veränderung der Geschwindigkeit), Dichte ρ, Viskosität γ, Aortenradius a, Radius der Klappenöffnung r und Klappenflügellänge l. Diese Werte lieferten die Differentialgleichung

p_{1} − p_{2} = −\frac{ρla}{r} · \frac{∂u}{∂t} − \frac{1}{2} · ρu^{2}(\frac{a^{4}}{r^{4}} − 1) .

Die berechneten Werte stimmten in erstaunlich hohem Maße mit den gemessenen überein.

Obwohl dem Konzept von Henderson und Johnson keine Schwachstellen nachgewiesen werden konnten, bevorzugten die meisten die ”Wirbel-Theorie“ von Bellhouse und Bellhouse, die der ”Breaking-of-a-jet-Theorie“ nicht grundlegend widerspricht. Es ist vielmehr eine andere Sichtweise bzw. ein anderes Modell für dieselben fluiddynamischen Gegebenheiten.

Ebenfalls im Jahre 1969 wagten B. J. Bellhouse und L. Talbot (siehe [4]) den ersten Versuch, die Funktion von Herzklappen mathematisch zu modellieren. Sie erreichten aber nur eine grobe Annäherung an die gemessenen Werte. Bessere Ergebnisse erzielte C. S. Peskin (siehe [30]) im Jahre 1977, indem er vollständige Bewegungsgleichungen aufstellte und diese analytisch löste. Um die Lösbarkeit zu gewährleisten, mußte er sich auf ein 2-dimensionales Modell beschränken und, um den Rechenaufwand in einem realisierbaren Rahmen zu halten, auf eine recht kleine Reynoldszahl (Maß für die Viskosität der betrachteten Flüssigkeit), die den physiologischen Gegebenheiten nicht gerecht wurde.

Die ersten akzeptablen Werte durch Lösen von fluiddynamischen Gleichungen erhielten C. S. F. Lee und L. Talbot 1979 (siehe [28] und [15]), indem sie die Gleichungen von 1969 noch weiter verfeinerten. Sie versuchten, den Effekt der Blutviskosität als Randbedingungen für den Prozeß der Klappenöffnung einzubringen, d.h. die Geschwindigkeit an verschiedenen Stellen der Klappenflügel sollte gleich sein. Zusammen mit den Gleichgewichtsgleichungen für Masse, Impuls und Energie erhielten sie eine übersichtliche aber doch genaue Beschreibung des Öffnugsprozesses. Für das Schließen der Klappen wurden die Gleichungen von 1969 übernommen, da dieser Vorgang nur unwesentlich durch Viskosität beeinflußt wird.

3.3 Ein moderner Ansatz: Doktorarbeit von Mary J. King, 1994

Im Jahre 1994 veröffentlichte Mary J. King an der Universität Leeds ihre Arbeit mit dem Titel ”Numerische und experimentelle Studien zum Fluß durch zweiflügelige mechanische Herzklappen“ und bewirkte damit einen großen Schritt in der detaillierten Berchnung der Herzklappenfunktion. Bis zu diesem Zeitpunkt existierte nur ein 2-dimensionales, laminares, zeitunabhängiges Modell bezüglich einer physiologisch nicht angemessenen Reynoldszahl, das den Verlauf der Stromlinien zwar grob annähern konnte, aber nicht in der Lage war, die Bildung von Wirbeln und Turbulenzen exakt vorher zu sagen. Das Ziel der Doktorarbeit von King war also ein 3-dimensionales, zeitabhängiges Modell einer zweiflügeligen Klappe in aortischer Position bezüglich einer physiologisch akzeptablen Reynoldszahl.

Als Grundlage, bzw. als experimentelle Vergleichswerte, wurden die mechanischen Daten der Carbomedics und der St. Jude Doppelflügelklappen benutzt und mit zwei leicht veränderten Prototypen verglichen. Der wesentliche Unterschied der vier Klappentypen liegt in den verschiedenen Öffnungswinkeln der Klappenflügel zwischen 78^{∘} und 85^{∘}.

Wie in der Doktorarbeit (vergl. [23], S. 6ff) ausführlich beschrieben, wird zum vereinfachten Verständnis der fluiddynamischen Verhältnisse zunächst ein Fluß um einen 2-D Zylinder herum angenommen. Bei einer sehr niederigen Reynolds-Zahl bewegt sich das Fluid wie in Abb. 11, A, um den Zylinder herum. Wird die Reynolds-Zahl erhöht, lösen sich die Randschichten ab und es entstehen zwei zusammenhängende Wirbel, wie man in Abb. 11, B, sehen kann. Diese Wirbel dehnen sich aus und werden instabil, sobald die Reynolds-Zahl noch weiter ansteigt (siehe Abb. 11, C), bis sie schließlich alternierend von gegenüberliegenden Seiten des Zyliders strömen. Werden Winkelgeschwindigkeit und Scherkräfte der Wirbel zu groß, so brechen sie zusammen und werden zu Turbulenzen (siehe Abb. 11, D). Als nächstes wird ei-

Abbildung 11: Entwicklung des Flußverhaltens um einen 2-D zylinder bei wachsender Reynolds-Zahl (aus [23], S.7)

ne Doppelflügel-Klappenprothese betrachtet. Die Klappenflügel bestehen aus semielliptischen flachen Platten, die, wenn die Klappe vollständig geöffnet ist, in einem bestimmten Winkel (siehe Abb. 12) zur Strömung geneigt sind. Da die Klappenflügel sich über die Klappenöffnung ausdehnen, kann man weitere Modelle auf eine zweidimensionale Platte in der Mitte der Flügel beschränken. Diese dünne Platte würde sich innerhalb des Flüssigkeitsstromes ähnlich verhalten wie der oben beschriebene Zylinder.

Der Herz-Zyklus wird des weiteren mit folgenden Annäherungs- bzw. Durchschnittswerten beschrieben:

• Die Flußrate Q wird durch eine Sinus-Kurve Q = Q_{m} sin ωt approximiert, wobei Q_{m} die maximale Flußrate ist. Durch ω können unterschiedliche Herzfrequenzen eingehen.

• Das ausgeworfene Blutvolumen pro Zyklus wird auf einen Durchschnittswert von 70 ml bei einer Herzfrequenz von 72 bpm (Schlägen pro Minute) festgelegt. Integriert man über den Abschnitt einer Herz-Zyklus-Kurve von 0 bis π, so läßt sich die maximale Flußrate und somit auch die Flußgeschwindigkeit berechnen.

• Man nimmt eine dynamische Blutviskosität von η = 4 · 10^{−3} P a · s, die Dichte von Blut mit ρ = 1000 kg · m^{−3}, eine menschliche Durchschnitts-Aorta (achsensymmetrisch bzw. gerades Rohr) mit r = 0, 0145 m Radius und eine Durchschnittsgeschwindigkeit des Blutes in der Aorta von u = 0, 3995 m · s^{−1} an und erhält so bzgl. des Durchmessers eine Reynoldszahl von Re_{D} = 3000, bzw. Re_{R} = 1500 bzgl. des Radius.

Re_{R} = \frac{ρur}{η} = 1000 · 0, 3995 · 0, 0145 · 250 ≈ 1500 .

Für Rohrströmung liegt diese Reynoldszahl zwischen den Werten für laminaren und turbulenten Fluß. Daher würde man für das betrachtete 2D Herzklappenmodell ein Flußverhalten erwarten, das zwischen den Ansätzen ’C’ und ’D’ in Abb. 11 liegt. Die genaue Bedeutung und Berechnung der Reynoldszahl Re wird in Abschnitt 4.3 noch detailliert beschrieben. Obwohl man Blut eigentlich nicht als Newtonsches Fluid bezeichnen kann, wird zur weiteren Modellierung von laminarer Strömung ausgegangen. Da sich das Blut nur im Zentrum der großen Arterien nicht wie ein Newtonsches Fluid verhält, ist diese Vereinfachung im Bereich direkt hinter einer Herzklappe, wo ungewöhnlich große Schubkräfte herrschen, durchaus gerechtfertigt.

Eine flache, leicht geneigt Fläche, wie in diesem Modell, wird sich im Blutstrom wie eine Flugzeugtragfläche im Luftstrom verhalten, d.h. sie wird Auftrieb erzeugen. Zur Modellierung von Doppelflügelprothesen können die Flügel also wie Tragflächen behandelt werden, die sich in einem Winkel zwischen 5^{∘} und 12^{∘} bewegen, wie in Abb. 12 angedeutet.

Betrachtet man ein ideales Fluid im Fluß über eine Tragfläche, erkennt man, daß der Oberflächendruck an der Unterseite der Tragfläche größer ist als der an der Oberseite, so daß eine Auftriebskraft entsteht. Dieser Druckunterschied kann z.B. auf eine Zirkulation um den Flügel herum zurückgeführt werden.

Die Theorie zu solchen Problemen liefert nun die Vorhersage von Wirbelbildungen an der Tragfläche, was auf die Klappenflügel übertragbar sein sollte. Um diese Wirbelbildungen genau erfassen zu können, muß besonderer Augenmerk auf die Viskosität des Fluids und somit auch auf die Trägheit der Grenzschichten gerichtet werden, die sich in der Flußgeschwindigkeit bemerkbar macht: Im Haupteil des Fluids findet man eine Geschwindigkeit u ≠ 0, während man an der Oberfläche der Tragfläche und an der Außenwand des Rohres aufgrund von Reibung als Randbedingung eine Geschwindigkeit von u = 0 voraussetzt.

Abbildung 12: schematische Darstellung einer Doppelflügelklappe, Modell vergleichbar mit einer Tragfläche (aus [23], S.8)

Aus diesen Voraussetzungen erhält man eine speziell angepaßte Variante der Navier-Stokes-Gleichungen mit dazugehörigen Randbedingungen. Die Software, die zur Lösung benutzt wurde, basiert auf der Methode finiter Elemente, benutzt gewichtete Residuen und die Galerkin-Methode.

Das heißt, der zu betrachtende Raum wird zerteilt, z.B. trianguliert, die entstandenen Knoten werden bezeichnet und das Problem mit den dazugehörigen Randbedingungen für die neuen Elemente spezialisiert. Die daraus resultierende globale Matrix bzw. das Gleichungssystem wurde iterativ gelöst, wobei bei den Toleranzgrenzen zwischen linearen und nicht-linearen Lösungen unterschieden wurde.

Die rein mathematischen und numerischen Fragestellungen, die diese Problematik der Navier-Stokes-Gleichungen hervorruft, sollen in den folgenden Abschnitten detailliert erläutert werden.

4 Navier-Stokes-Gleichungen und Reynoldszahl

Um die Navier-Stokes-Gleichungen herleiten und verstehen zu können (vergl. [7] und [38]), macht man häufig einen Umweg und betrachtet zuerst ein ideales Fluid, d. h. eine Flüssigkeit ohne innere Reibung bzw. ohne Zähigkeit, und ergänzt anschließend in der beschreibenden Gleichung Ausdrücke, die die Reibung beschreiben. Zur genauen Markierung von verschieden Orten innerhalb des Fluids zu verschiedenen Zeiten versieht man das Modell mit einem Koordinatensystem, in dem die Bahnlinie einzelner Punkte in Raum und Zeit verfolgt werden kann.

Im weiteren werden folgende Bezeichnungen verwendet:

• D sei ein Gebiet in ℝ^{2} oder ℝ^{3} gefüllt mit Flüssigkeit.

• \underline{x} ∈ D sei ein Punkt in D; betrachtet werden die Flüssigkeits-Teilchen, die zum Zeitpunkt t den Punkt \underline{x} passieren.

• \underline{u}(\underline{x}, t) sei die Geschwindigkeit der Teilchen, die \underline{x} zum Zeitpunkt t passieren.

• Man nimmt an, daß das Fluid eine zu jedem Zeitpunkt t und an jedem Ort \underline{x} wohldefinierte Dichte ρ(\underline{x}, t) > 0 hat.

• p(\underline{x}, t) beschreibt den orts- und zeitabhängigen Druck, \underline{b}(\underline{x}, t) die auf einen Punkt wirkende Schwerkraft.

• W ⊂ D Teilmenge mit glattem Rand, so ist die Masse m(W, t) der Flüssigkeit in W zum Zeitpunkt t gegeben durch

m(W, t) = ∫_{W} ρ(\underline{x}, t) dV ,

wobei dV ein Volumen- bzw. Flächenelement beschreibt.

4.1 Herleitung der Euler-Gleichungen

Die Herleitung der Euler-Gleichungen basiert auf drei elementaren, physikalischen Gesetzen:

1. Massenerhaltung

2. Impulserhaltung

3. Energieerhaltung

Zu 1.: In einem beliebig abgegrenzten Volumen W ⊂ D, welches fest gewählt und unabhängig von der Zeit sei, muß die zeitliche Abnahme der darin enthaltenen Masse gleich der Differenz von ein- und ausströmender Masse sein. Diese Massenveränderung wird beschrieben durch:

\frac{∂}{∂t}m(W, t) = \frac{∂}{∂t} ∫_{W} ρ(\underline{x}, t) dV = ∫_{W} \frac{∂ρ}{∂t}(\underline{x}, t) dV .

Die Differenz des Massenflusses durch ein Oberflächenelement dA des Volumens W in Richtung der äußeren Normalen \underline{n} von A erhält man über das Oberflächenintegral ∮ ρ\underline{u} · \underline{n} dA. Das Gesetzt der Massenerhaltung lautet also

∫_{W} \frac{∂ρ}{∂t}(\underline{x}, t) dV = − ∮_{∂W} ρ\underline{u} · \underline{n} dA .

Die Anwendung des Gaußschen Integralsatzes liefert

∫_{W} \frac{∂ρ}{∂t}(\underline{x}, t) dV = − ∮_{∂W} ρ\underline{u} · \underline{n} dA =^{Gauß} − ∫_{W} div(ρ\underline{u}) dV .

Da W ⊂ D beliebig gewählt war, führt dies zu

∫_{W} \{\frac{∂ρ}{∂t} + div(ρ\underline{u})\} dV = 0 ⟺ \frac{∂ρ}{∂t} + div(ρ\underline{u}) = 0 .  (1)

(Differentialform der Massenerhaltung bzw. Kontinuitätsgleichung)

Zu 2.: Man unterscheidet stets zwischen zwei Arten von Kräften, die auf ein fest gewähltes Volumenteil W des Fluids wirken können:

Einerseits Druck auf die Oberfläche, also direkte Krafteinwirkung durch benachbarte Materie, und andererseits externe Kräfte wie Schwerkraft oder durch magnetische Felder verursachte Kräfte. Da wir zunächst nur eine ideale Flüssigkeit, d. h. eine Flüssigkeit ohne innere Reibung, betrachten, findet man an der Oberfläche ∂W keine Tangentialkräfte, sondern ausschließlich orthogonale Kräfte in Richtung der äußeren Normalen. Wie unter 1. ist der Druck p bzw. die Kraft \underline{S}_{∂W}, die in Richtung eines äußeren Normalenvektors auf den Rand von W wirkt, mittels Gaußschem Integralsatz gegeben durch

\underline{S}_{∂W} = − ∮_{∂W} p · \underline{n} dA = − ∫_{W} gradp dV .  (2)

Wenn \underline{b}(\underline{x}, t) die gegebene Anziehungskraft pro Masseneinheit bezeichnet, so ist analog die gesamte Schwerkraft

\underline{B} = ∫_{W} ρ\underline{b} dV ,

so daß insgesamt auf jede Volumeneinheit folgende Kräfte einwirken:

Kraft pro Volumeneinheit = −gradp + ρ\underline{b} .

Bezeichnet man mit \underline{x}(t) = (x(t), y(t), z(t)) die Kurve, die ein Fluidpartikel beschreibt, so läßt sich dessen Geschwindigkeit \underline{u} schreiben als \underline{u}(\underline{x}(t), t) = \frac{d\underline{x}}{dt}(t), so daß man für die Beschleunigung folgende Schreibweise erhält:

\frac{d\underline{u}}{dt}(\underline{x}(t) , t) = \frac{d}{dt}\underline{u}(x(t), y(t), z(t), t)
= \frac{∂\underline{u}}{∂x} · ẋ + \frac{∂\underline{u}}{∂y} · ẏ + \frac{∂\underline{u}}{∂z} · ż + \frac{∂\underline{u}}{∂t}
=: grad\underline{u} · \underline{u} + ∂_{t}\underline{u}
=: \frac{D\underline{u}}{Dt} , wobei \frac{D}{Dt} := \underline{u} · \underline{∇} + ∂_{t} .  (3)

Mit dem 2. Gestz von Newton (Kraft = Masse × Beschleunigung) gelangt man von hier zur Differentialform der zweiten Gleichung:

ρ · \frac{D\underline{u}}{Dt} = −gradp + ρ\underline{b}  (4)

(Impulserhaltung)

Ein Fluid heißt inkompressibel, wenn für eine beliebige Teilmenge W des Fluids gilt, daß das Volumen ∫_{W} dV von W konstant in t ist, woraus man die Bedingung div \underline{u} = 0 herleiten kann, falls das Fluid von einer festen Wand umgeben ist.

Aus der Kontinuitätsgleichung

\frac{∂ρ}{∂t} + div(ρ\underline{u}) = ∂_{t}ρ + gradρ · \underline{u} + ρdiv\underline{u}
=: \frac{Dρ}{Dt} + ρdiv\underline{u}
= 0  (5)

und der Tatsache, daß stets ρ > 0 ist, folgt, daß die Kontinuitätsgleichung für ein inkompressibles Fluid

\frac{Dρ}{Dt} = 0  (6)

lautet.

Man kann also sagen, daß ein Fluid inkompressibel ist, wenn die Massen-Dichte ρ konstant den Bewegungen des Fluids folgt. Ist ρ nicht nur unabhängig von der Zeit, sondern zusätzlich noch unabhängig vom Ort, so spricht man von einem homogenen Fluid. Eine solche Flüssigkeit, die zum Zeitpunkt t = 0 homogen ist, muß nicht zwingend homogen bleiben; dieses ist nur dann der Fall, wenn sie auch inkompressibel ist.

Zu 3.: Betrachtet man den Flüssigkeitsstrom in drei Dimensionen, so erhält man mit \underline{u}(\underline{x}(t), t) = \frac{d\underline{x}}{dt}(t) = (ẋ, ẏ, ż) (t) aus den bisherigen Gleichungen ein Gleichungssystem mit 4 Gleichungen, aber 5 Unbekannten:

1. \frac{Dρ}{Dt} + ρdiv\underline{u} = 0 (Masse)

2. ρ · \frac{D\underline{u}}{Dt} = −gradp + ρ\underline{b} (Impuls)

Die noch fehlende Gleichung gewinnt man aus dem Gesetz der Energieerhaltung.

Man nimmt an, daß die Gesamtenergie des Fluids geschrieben werden kann als:

E_{total} = E_{kin} + E_{int} + E_{ext} ,

wobei die kinetische Energie in einem Volumen D gegeben ist durch

E_{kin} = \frac{1}{2} ∫_{D} ρ∥\underline{u}∥^{2} dV .

E_{int} bezeichnet die innere Energie bzw. die Wärme, die sich aus intermolekularen Potentialen, molekularen Schwingungen u. ä. zusammensetzt. Die aus externen Kräften, wie etwa der Schwerkraft, resultierende Energie E_{ext} = \underline{b} wird meist vernachlässigt, da sie nur minimal ins Gewicht fällt. Aus diesen Voraussetzungen kann man verschiedene Varianten von Energiegleichungen ableiten. Nimmt man z. B. E_{int} als konstant an, so muß E_{kin} unabhängig von Bewegung bzw. Zeit sein, um das Gleichgewicht zu erhalten, d. h. es handelt sich um einen inkompressibeln Fluß. Es gilt also

\frac{∂}{∂t} (\frac{1}{2} ∫_{D} ρ∥\underline{u}∥^{2} dV ) = 0 .  (7)

Für jede stetig differenzierbare Funktion f(\underline{x}, t) gilt

\frac{∂}{∂t} ∫_{D} f dV = ∫_{D} \frac{Df}{Dt} dV .

Dies liefert zusammen mit der Kettenregel und der Kontinuitätsgleichung für inkompressiblen Fluß (6):

0 = \frac{∂}{∂t} (\frac{1}{2} ∫_{D} ρ∥\underline{u}∥^{2} dV )

= \frac{1}{2} ∫_{D} (\underbrace{\frac{Dρ}{Dt}}_{=0} ·∥\underline{u}∥^{2} + 2ρ\underline{u}\frac{D\underline{u}}{Dt}) dV

= ∫_{D} ρ · \underline{u} · \frac{D\underline{u}}{Dt} dV = 0 .  (8)

Das Gleichgewicht des Momentums liefert mit \underline{b} = 0 weiter

− ∫_{D} \underline{u} · gradp dV = 0 .  (9)

Nimmt man sinnvollerweise an, daß das Fluid den äußeren Rand von D nicht überschreitet, also \underline{u} · \underline{n} = 0 auf ∂D, und integriert partiell, so folgt

− ∫_{D} \underline{u} · gradp dV = \underbrace{− ∮_{∂D} \underline{u} · p · \underline{n} dA}_{=0} + ∫_{D} div(\underline{u})p dV

= ∫_{D} div(\underline{u})p dV = 0  (10)
⟹ div \underline{u} = 0 ,  (11)

da der Fall p = 0, d.h. nicht vorhandener Druck, für die hier betrachtete Fragestellung physikalisch nicht sinnvoll ist.

Also haben die Euler-Gleichungen für ein inkompressibles Fluid insgesamt folgende Gestalt:

\begin{matrix}ρ · \frac{D\underline{u}}{Dt} & = & −gradp & (Impuls) \\ \frac{Dρ}{Dt} & = & 0 & (Masse) \\ div \underline{u} & = & 0 & (Energie) \end{matrix} \} auf D

mit \underline{u} · \underline{n} = 0 auf ∂D .  (12)

4.1.1 Beispiel

Schon an einem einfachen Beispiel kann man erkennen, daß die Euler-Gleichungen wegen der fehlenden Reibungskomponente nicht in der Lage sind, das Flußverhalten eines realen Fluids vollständig zu beschreiben.

Betrachtet man ein inkompressibles, homogenes Fluid beim Fluß durch eine ”zweidimensionale“ Röhre in x-Richtung, so ergibt sich ein deutlicher Widerspruch. Gegebene Voraussetzungen sind also:

• ρ = ρ_{0} konstant, da das Fluid inkompressibel und homogen ist; die Kontinuitätsgleichung entfällt also.

• Da eine Lösung der Form

\underline{u}(x, y, t) = (\underline{u}(x, t), 0)
und p(x, y, t) = p(x)

gesucht ist, liefert die Energie-Gleichung

div\underline{u} = 0 ⇒ ∂_{x}\underline{u} = 0
⇒ \underline{u} · grad\underline{u} = 0 .

• Dieses vereinfacht auch die Gleichung zur Impulserhaltung zu

ρ_{0}(∂_{t}\underline{u} + \underline{u} · grad\underline{u}) = ρ_{0} · ∂_{t}\underline{u} = −grad p = ∂_{x}p .  (13)

Man beschränkt sich nun auf ein Rohstück in x-Richtung von 0 − L und bezeichnet die Druckverhältnisse am Rand mit p(0) = p_{1} und p(L) = p_{2}, wobei natürlich p_{1} > p_{2} sein muß.

[figure: Flußrichtung; Länge L; \begin{matrix}Druck \\ = p_{1} \end{matrix}; \begin{matrix}Druck \\ = p_{2} \end{matrix}]

Für p erhalten wir also aus Gleichung (13):

ρ_{0}∂_{t}\underline{u} = −∂_{x}p
⇒ ∂_{x}^{2}p = 0 mit p(0) = p_{1} , p(L) = p_{2}

⇒ p(x) = p_{1} − (\frac{p_{1} − p_{2}}{L} ) · x
⇒ −∂_{x}p = \frac{p_{1} − p_{2}}{L} = ρ_{0}∂_{t}\underline{u}
⇒ \underline{u} = \frac{p_{1} − p_{2}}{ρ_{0} · L} · t + konst.

D. h. bei konstantem Druckgradienten ∂_{x}p bzw. konstantem Druckunterschied p_{1}−p_{2} wächst die Geschwindigkeit \underline{u} unendlich an. Dieser Widerspruch kann nur entstehen, weil die Reibung nicht berücksichtigt wurde.

4.2 Herleitung der Navier-Stokes-Gleichungen

Im folgenden werden einige Ausdrücke zur besseren Übersichtlichkeit mit Nablabzw. Laplace-Operatoren bezeichnet:

Den Laplace-Operator Δ eines skalaren Feldes u(D) erklärt man für kartesische Koordinaten durch den skalaren Wert

Δ u = \frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} + \frac{∂^{2}u}{∂z^{2}}

Ist v(D) = (v_{x}, v_{y}, v_{z})(D) ein Vektorfeld, so definiert man ebenfalls für kartesische Koordinaten Δv durch den Vektor

Δ v = [(Δ v_{x}), (Δ v_{y}), (Δ v_{z})]^{t}

Der Nabla-Operator ∇ = (\frac{∂}{∂x}, \frac{∂}{∂y}, \frac{∂}{∂z}) hingegen bezeichnet für ein skalares Feld u(D) den Vektor ∇ u = grad u und für ein Vektorfeld v(D) = (v_{x}, v_{y}, v_{z})(D) andererseits entweder den skalaren Wert ∇ v = div v, oder die Jacobi-Matrix von \underline{u} mit der Spur div \underline{u}.

Betrachtet man eine Ebene S innerhalb des Fluids, so wurde bis jetzt vorausgesetzt, daß alle Kräfte die auf S wirken, orthogonal zu dieser sind. Dieser Gedanke sollte nun um Tangential- bzw. Reibungskräfte erweitert werden.

Statt wie bisher anzunehmen, daß gilt

Kraft auf S pro Flächeneinheit = p(\underline{x}, t) · \underline{n}

verändert man diese Definition zu

Kraft auf S pro Flächeneinheit = p(\underline{x}, t) · \underline{n} + \underline{\underline{σ}}(\underline{x}, t) · \underline{n} ,

wobei die Reibungskomponente oder der Tensor \underline{\underline{σ}} eine Matrix mit 4 wichtigen Eigenschaften ist:

1. \underline{\underline{σ}} · \underline{n} muß nicht zwingend parallel zu \underline{n} sein.

2. \underline{\underline{σ}} hängt ausschließlich von der Jacobi-Matrix ∇\underline{u} ab und ist proportional zu diesem Gradienten.

3. \underline{\underline{σ}} ist symmetrisch.

4. \underline{\underline{σ}} ist invariant unter Körper-Rotation und Translation, d. h. für jede orthogonale Matrix \underline{\underline{U}} gilt

\underline{\underline{σ}}(\underline{\underline{U}} · ∇\underline{u} · \underline{\underline{U}}^{−1}) = \underline{\underline{U}} · \underline{\underline{σ}}(∇\underline{u}) · \underline{\underline{U}}^{−1} .

Aus 2. und 3. folgt:

\underline{\underline{σ}} hängt nur von dem symmetrischen Teil von ∇\underline{u} ab. Mit ∇\underline{u} = \underline{\underline{D}} + \underline{\underline{S}} und

\underbrace{\underline{\underline{D}} =: \frac{1}{2}[∇\underline{u} + (∇\underline{u})^{t}]}_{Diagonale} , \underbrace{\underline{\underline{S}} =: \frac{1}{2}[∇\underline{u} − (∇\underline{u})^{t}]}_{Rest}

folgert man nun weiter:

\underline{\underline{σ}} ist linear abhängig von \underline{\underline{D}} .
⇒ \underline{\underline{σ}} und \underline{\underline{D}} kommutieren, d.h. \underline{\underline{σ}} · \underline{\underline{D}} = \underline{\underline{D}} · \underline{\underline{σ}} .
⇒ \underline{\underline{σ}} und \underline{\underline{D}} können simultan diagonalisiert werden.
⇒ Die Eigenwerte von \underline{\underline{σ}} sind Linearkombinationen der Eigenwerte von \underline{\underline{D}} .
⇒^{4.} Die Eigenwerte von \underline{\underline{σ}} müssen ”symmetrisch“ sein, da man \underline{\underline{U}}
so wählen kann, daß zwei Eigenwerte von \underline{\underline{D}} permutiert werden.

Die einzigen in diesem Sinne symmetrischen Linearkombinationen sind

σ_{i} = λ(d_{1} + d_{2} + d_{3}) + 2μd_{i} mit i = 1, 2, 3 ,

wobei σ_{i} die Eigenwerte von \underline{\underline{σ}} und d_{i} die Eigenwerte von \underline{\underline{D}} bezeichnen. Der Parameter μ wird als 1. Viskositätskoeffizient bezeichnet.

Da die Spur der Jacobi-Matrix von \underline{u} = (u_{1}, u_{2}, u_{3}) gegeben ist durch

d_{1} + d_{2} + d_{3} = ∂_{x}u_{1} + ∂_{y}u_{2} + ∂_{z}u_{3} = div\underline{u}

läßt sich \underline{\underline{σ}} also als Linearkombination aus \underline{\underline{D}} und (div \underline{u})\underline{\underline{I}} schreiben, wobei \underline{\underline{I}} die Einheitsmatrix bezeichnet.

\underline{\underline{σ}} = λ(div\underline{u})\underline{\underline{I}} + 2μ\underline{\underline{D}}
=: 2μ[ \underline{\underline{D}} − \frac{1}{3}(div\underline{u})\underline{\underline{I}}] + ξ(div\underline{u})\underline{\underline{I}}  (14)

mit μ ˆ= 1. Viskositätskoeffizient
und ξ = (λ + \frac{2}{3}μ) ˆ= 2. Viskositätskoeffizient

Aus der ursprünglichen Gleichung zur Impulserhaltung folgt also mit dem 2. Gesetz von Newton (Masse × Beschleunigung = Kraft)

∫_{W} ρ · \frac{D\underline{u}}{Dt} dV = ∫_{∂W} (−p · \underline{n} + \underline{\underline{σ}} · \underline{n}) dA
(Gauß) = ∫_{W} (−∇p + ∇\underline{\underline{σ}}) dV
wobei ∇\underline{\underline{σ}} = ∇(2μ[ \underline{\underline{D}} − \frac{1}{3}(div\underline{u})\underline{\underline{I}}] + ξ(div\underline{u})\underline{\underline{I}})
= 2μ[ ∇\underline{\underline{D}} − \frac{1}{3}∇(div\underline{u}) ] + ξ∇(div\underline{u})
= 2μ(\frac{1}{2}Δ\underline{u} + \frac{1}{2}∇(div\underline{u}) − \frac{1}{3}∇(div\underline{u})) + ξ∇(div\underline{u})
= μΔ\underline{u} + (λ + μ)∇(div\underline{u})

⇒ ∫_{W} ρ\frac{D\underline{u}}{Dt} dV = ∫_{W} (−∇p + μΔ\underline{u} + (λ + μ)∇(div\underline{u})) dV
⇒ ρ\frac{D\underline{u}}{Dt} = −∇p + μΔ\underline{u} + (λ + μ)∇(div\underline{u}) .  (15)

Zusammen mit der Kontinuitätsgleichung und einer Ernergieerhaltungsgleichung beschreibt Gleichung (15) vollständig die Flußbewegung eines viskosen Fluids. Im Falle eines inkompressiblen, homogenen Flusses mit ρ ≡ ρ_{0} = konst. reduziert sich das komplette Gleichungssystem auf folgendes:

\frac{D\underline{u}}{Dt} = −gradp̃ + νΔ\underline{u}
mit p̃ = p/ρ_{0} , ν = μ/ρ_{0}
div \underline{u} = 0  (16)

Navier-Stokes-Gleichungen für inkompressiblen, homogenen Fluß

Diese Gleichungen müssen noch durch Randbedingungen ergänzt werden. Die passende Randbedingung zu den Euler-Gleichungen war \underline{u} · \underline{n} = 0, d.h. der Fluß kann den Rand nicht kreuzen, sich aber tangential dazu bewegen. Die Navier-Stokes-Gleichungen enthalten nun den zusätzlichen Term νΔ\underline{u}, durch den neben den 1. Ableitungen von \underline{u} nun auch die 2. Ableitungen eingebracht werden. Sowohl aus experimentellen, als auch aus mathematischen Gründen fordert diese Tatsache direkt eine 2. Randbedingung, da das Ziel eine eindeutige, stetig von den Anfangswerten abhängende Lösung ist.

Beobachtet man die Grenzschichten einer fließenden viskosen Flüssigkeit in einem Rohr, so kann man erkennen, daß auch die Tangentialgeschwindigkeit gegen 0 geht, was also insgesamt die Randbedingung

\underline{u} = 0 auf ∂D ,

natürlich nur bei fester, unbeweglicher Wand, ergibt.

4.3 Die Reynolds-Zahl

Die oben hergeleiteten Gleichungen reichen zwar zum vollständigen Beschreiben aller Flüssigkeitsbewegungen aus, sind jedoch häufig nicht analytisch zu lösen. Die wenigen berechenbaren Strömungen haben die Eigenschaft, daß alle Strömungsgrößen wie Geschwindigkeit, Druck usw. in jedem Raum-Zeit-Punkt des Strömungsfeldes definiert und berechenbar sind. Bahn- und Stromlinien bilden glatte Raumkurven, d. h. die Strömung besteht in einem Vorübergleiten benachbarter Flüssigkeitsschichten. Man spricht deshalb hier von laminarer Strömung.

In den meisten für die Praxis interessanten Fällen erweist sich diese Laminarströmung aber als instabil; es ensteht eine völlig andersartige turbulente Strömungsform.

Hier sind dann selbst bei unveränderlichen, stationären Randbedingungen Größen wie Geschwindigkeit und Druck im festgehaltenen Raumpunkt nicht mehr zeitlich konstant, sondern schwanken schnell und völlig unregelmäßig um einen zeitlichen Mittelwert, der sich bei instationärer turbulenter Strömung sogar noch ändern kann.

Da die Gleichung zur Impulserhaltung für turbulente Strömung aus den obigen Gründen differiert, ist es notwendig, durch einfache Modellversuche festzustellen, ob es sich in dem zu untersuchenden Fall um eine laminare oder turbulente Strömung handelt. Ein Maß für die Turbulenz ist die Reynolds-Zahl Re, die den Effekt von Viskosität auf das Flußverhalten beschreibt.

Für ein gegebenes Problem werden folgende charakterisierenden Größen eingeführt:

L ˆ= charakteristische Länge
U ˆ= charakteristische Geschwindigkeit
T = L/U ˆ= betrachtete Zeit

Dies ist eine äußerst ungenaue Definition, da z. B. bei einem Fluß durch einen Zylinder sowohl der Radius, als auch der Durchmesser als charakteristische Länge betrachtet werden könnte. Mittels dieser für das Problem typischen Größen führt man nun dimensionslose oder auch ”normierte“ Variablen ein:

˜\underline{u} := \underline{u}/U , ˜ \underline{x} := \underline{x}/L , t̃:= t/T

Die x-Komponente der Navier-Stokes-Gleichung für inkompressiblen Fluß lautet für \underline{u} = (u, v, w) und ˜ \underline{u} = (ũ, ṽ, w̃)

\frac{∂ u}{∂t} + u\frac{∂ u}{∂x} + v \frac{∂ u}{∂y} + w\frac{∂ u}{∂z} = −\frac{1}{ρ_{0}} \frac{∂p}{∂x} + ν [\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} + \frac{∂^{2}u}{∂z^{2}}] .  (17)

Ein Tausch der Variablen führt zu

\frac{∂ũ}{∂t̃} + ũ\frac{∂ũ}{∂x̃} + ṽ\frac{∂ũ}{∂ỹ} + w̃\frac{∂ũ}{∂z̃} = −\frac{1}{ρ_{0}U^{2}} \frac{∂p}{∂x̃} + \frac{ν}{L · U} [\frac{∂^{2}ũ}{∂x̃^{2}} + \frac{∂^{2}ũ}{∂ỹ^{2}} + \frac{∂^{2}ũ}{∂z̃^{2}}] .  (18)

Mit p̃:= p/ρ_{0}U^{2}, der kinetischen Viskosität ν := μ/ρ_{0} und den analogen y- und z-Komponenten haben die dimensionslosen Navier-Stokes-Gleichungen folgende Gestalt:

\frac{∂ ˜ \underline{u}}{∂ t̃} + (˜\underline{u} · ∇) ˜ \underline{u} = −grad p̃ + \frac{μ}{ρ_{0} L U} Δ ˜ \underline{u}
div ˜ \underline{u} = 0 auf D
˜\underline{u} = 0 auf ∂D .  (19)

Man definiert die Reynolds-Zahl Re nun durch

Re := \frac{charakteristische Geschwindigkeit × Länge}{kinematische Viskosität}
= \frac{L · U}{ν}
= \frac{ρ_{0} · L · U}{μ} .

Sie kennzeichnet die Größenordnung des Verhältnisses von Trägheits- zu Zähigkeitskräften. Eine kleine Re-Zahl entspricht geringer Geschwindigkeit, sehr zäher Flüssigkeit oder kleinen räumlichen Ausmaßen und steht damit für ein Überwiegen der Zähigkeitskräfte. Bei großer Re-Zahl spielt die Zähigkeit eine kleinere Rolle und es kommt vor allem auf das Gleichgewicht zwischen Trägheits-, Druck- und äußeren Kräften an. Die in der Praxis interessanteren Fälle sind also diejenigen, die durch eine hohe Reynolds-Zahl ausgezeichnet sind.

Setzt man die Kraftdichte bzw. die externen Kräfte nicht vereinfacht auf 0, so müssen die obigen Gleichungen (19) noch um einen Summanden f = f(x, t) ergänzt werden, der eben diese Kräfte beschreibt. Auf diese Weise erhält man die allgemeinen Navier-Stokes-Gleichungen, wie sie in den meisten Standardwerken zur Strömungslehre zu finden ist:

\frac{∂ \underline{u}}{∂ t} + (\underline{u} · ∇) \underline{u} = −grad p + \frac{1}{Re} Δ \underline{u} + f
div \underline{u} = 0 auf D
\underline{u} = 0 auf ∂D

(20)

5 Theoretische Vorbereitungen

Differentialgleichungen können auf verschieden Arten numerisch oder analytisch behandelt bzw. gelöst werden (vergl. [36],[11], [17], [8], [31], [9] und [10]). Eine häufig benutzte Technik sind die Differenzenverfahren, deren Grundidee darin besteht, auftretende Ableitungen der gesuchten Funktion durch Informationen an diskreten Stellen mittels Gitterfunktionen zu approximieren, also Ableitungen durch geeignete Differenzenquotienten näherungsweise zu ersetzen. Die Methode, die im folgenden näher betrachtet werden soll, ist die Methode der finiten Elemente, die sich von den Differenzenverfahren grundlegend unterscheidet:

Anstelle der punktweisen Gültigkeit wird eine unter gewissen Zusatzbedingungen äquivalente Integral- oder Variationsbedingung zugrunde gelegt, die durch die Aproximationen u_{n} erfüllt werden muß. Bei diesen Approximationen handelt es sich nun nicht mehr um Gitterfunktionen, sondern um Summen von ”echten“ Funktionen, die nur begrenzt oft differenzierbar sind, dafür aber kleine Träger besitzen. Um die Lösbarkeit der Variationaufgabe zu sichern, müssen angepaßte Funktionenräume verwendet werden, die genauere Aussagen z.B. zum Definitionsbereich der Lösung und ihren Approximationen zulassen.

5.1 Funktionenräume

Die Navier-Stokes-Gleichungen machen nur dann mathematisch und physikalisch Sinn, wenn der betrachtete fluidgefüllte Raum Ω eine offene, aber beschränkte Teilmenge des ℝ^{n} ist mit n = 2 oder n = 3. Der Rand von Ω sei von nun an bezeichnet mit Γ, bzw. mit der Funktion Γ : ℝ^{n} → ℝ^{n−1}. Im weiteren wird immer wieder eine bestimmte ”Glattheitseigenschaft“ von Γ verlangt, die sich folgendermaßen darstellt (vgl.: [14], S. 102):

Definition 5.1 Ω ⊂ ℝ × ℝ^{n−1}. Sei Γ : Ω → ℝ^{n−1} eine Funktion. Man sagt, Γ genüge in ℝ^{n} lokal einer Lipschitz-Bedingung mit der Lipschitz-Konstanten L ≥ 0, falls jeder Punkt (x, y) ∈ Ω eine Umgebung U und eine davon abhängige Konstante L ≥ 0 besitzt, so daß für alle (x.y), (x, ỹ) ∈ Ω gilt:

∥Γ(x, y) − Γ(x, ỹ)∥ ≤ L∥y − ỹ∥ .

Im folgenden sei immer vorausgesetzt, daß Ω ⊂ ℝ^{n} offen und beschränkt mit Lipschitz-stetigem Rand ist.

Bezeichnungen Sei U ⊂ ℝ^{n} eine offene Menge. Wir bezeichnen mit C(U) den Vektorraum aller stetigen Funktionen f : U → ℝ und mit C_{0}(U) den Untervektorraum aller Funktionen f ∈ C(U), die kompakten Träger in U haben. Für eine natürliche Zahl k sei C^{k}(U) der Vektorraum der k-mal stetig partiell differenzierbaren Funktionen f : U → ℝ sowie

C^{∞}(U) := ⋂_{k=0}^{∞} C^{k}(U) ,

C_{0}^{k}(U) := C_{0}(U) ∩ C^{k}(U) .

Für die Beweise einiger Existenzsätze sind geeignete Hilberträume von Nöten, z.B. die Sobolevräume, deren Verständnis nur durch den Weg über die Lebesguesche Theorie gewährleistet werden kann. Im Unterschied zur Riemannschen Theorie werden Ober- und Unterintegrale hier mit Hilfe der halbstetigen Funktionen aus H^{↑}(ℝ^{n}) bzw. H^{↓}(ℝ^{n}) anstelle der Treppenfunktionen definiert.

Definition 5.2 Die Menge aller Funktionen f : ℝ^{n} → ℝ ∪ {∞}, die sich als Limiten monoton wachsender (bzw. fallender) Funktionenfolgen f_{ν} aus C_{0}(ℝ^{n}) darstellen lassen, d.h.

f(x) := lim_{ν→∞} f_{ν}(x) für alle x ∈ ℝ^{n} ,

bezeichnet man mit H^{↑}(ℝ^{n}) (bzw. H^{↓}(ℝ^{n})) oder kurz H^{↑} (bzw. H^{↓}).

Definition 5.3 (Oberintegral, Unterintegral) Sei f : ℝ^{n} → ℝ ∪ {±∞} eine beliebige Funktion. Dann setzt man

∫ ^{∗} f(x)dx := inf {∫_{ℝ^{n}} φ(x)dx : φ ∈ H^{↑}, φ ≥ f} ,

∫_{∗} f(x)dx := sup {∫_{ℝ^{n}} ψ(x)dx : ψ ∈ H^{↓}, ψ ≤ f} .

Wie in der Riemannschen Theorie heißt eine Funktion nun integrierbar, falls Oberund Unterintegral übereinstimmen.

Definition 5.4 Eine Funktion f : ℝ^{n} → ℝ ∪ {±∞} heißt Lebesgue-integrierbar, falls

−∞ < ∫_{∗} f(x)dx = ∫ ^{∗} f(x)dx < +∞ .

Der gemeinsame Wert des Ober- und Unterintegrals heißt dann das Lebesgue-Integral von f und wird mit ∫ f(x)dx bezeichnet.

Bezeichnung Wir bezeichnen mit L_{1}(ℝ^{n}) die Menge aller Lebesgue-integrierbaren Funktionen f : ℝ^{n} → ℝ, ohne die Werte ±∞. Mit den Rechenregeln der Integralrechnung wird L_{1}(ℝ^{n}) zu einem Vektorraum, auf dem das Integral ein lineares, monotones Funktional darstellt.

Der letzte Schritt auf dem Weg zu dem ersten Ziel, den Lebesgue-Räumen, ist die Definition meßbarer Mengen, bzw. meßbarer Funktionen:

Definition 5.5 Sei M ⊂ ℝ^{n} eine beliebige Teilmenge. Eine Funktion f : M → ℝ ∪ {±∞} heißt Lebesgue-meßbar auf M, falls die trivial fortgesetzte Funktion f̃ mit

f̃(x) := \{ \begin{matrix}f(x) für & x ∈ M \\ 0 für & x ∈ ℝ^{n}\M \end{matrix}

über ℝ^{n} Lebesgue-integrierbar ist. Man setzt dann

∫_{M} f(x)dx := ∫_{ℝ^{n}} f̃(x)dx .

Definition 5.6 (Lebesgue-Räume) Sei Ω eine offene Teilmenge des ℝ^{n}, n ∈ ℕ, p ∈ ℝ, m ≥ 0, p ≥ 1. Man bezeichnet mit L^{p}(Ω) den Lebesgueschen Raum. Die Elemente dieses Raumes sind Klassen äquivalenter Funktionen f(x), die auf Ω Lebesgue-meßbar sind und für die |f(x)|^{p} Lebesgue-summierbar ist, d.h.

∥f∥_{p} = [∫_{Ω} |f(x)|^{p}dx]^{1/p} ≠ 0 für 1 ≤ p < ∞
und ∥f∥_{∞} = inf{K : |f(x)| ≤ K fast überall} ≠ 0 .

Um letztendlich Sobolevräume definieren zu können, benötigen wir noch die Definition einer weiteren Norm. Mit α = (α_{1}, …, α_{n}) , α_{i} ∈ ℕ , |α| = α_{1} + … + α_{n} und

D^{α} = D_{1}^{α_{1}} ⋯ D_{n}^{α_{n}} = \frac{∂^{|α|}}{∂x_{1}^{α_{1}} ⋯ ∂x_{n}^{α_{n}}}

hat diese Norm folgende Gestalt:

∥u∥_{m,p} = [ ∑_{0≤|α|≤m}(∥D^{α}u∥_{p})^{p}]^{1/p} , 1 ≤ p < ∞

und ∥u∥_{m,∞} = max_{|α|≤m} ∥D^{α}u∥_{∞} .

Definition 5.7 (Sobolev-Räume durch Vervollständigung) Unter der Vervollständigung eines metrischen Raumes X versteht man einen über Cauchy-Folgen (wie bei der Vervollständigung von ℚ zu ℝ) gewonnenen vollständigen Raum X̃, der zu X isometrisch und isomorph ist. Auf diese Weise erhält man die erste Definition von Sobolev-Räumen.

H^{m,p}(Ω) := Vervollständigung von {u ∈ C^{∞}(Ω) : ∥u∥_{m,p} < ∞} bzgl. ∥ · ∥_{m,p}
H_{0}^{m,p}(Ω) := Vervollständigung von {u ∈ C_{0}^{∞}(Ω) : ∥u∥_{m,p} < ∞} bzgl. ∥ · ∥_{m,p}

Obwohl diese Definition von H^{m,p}(Ω) und H_{0}^{m,p}(Ω sehr einfach und häufig nützlich ist, hat sie den Nachteil, daß die funktionstheoretischen Eigenschaften der Elemente von H^{m,p}(Ω) nicht geklärt sind. Deshalb wird zusätzlich ein anderer Weg der Definition über schwache Ableitungen eingeschlagen.

Sei Ω offene Teilmenge des ℝ^{n}. Sei X := C_{0}^{∞}(Ω), also die Menge aller unendlich oft differenzierbaren Funktionen Φ auf Ω mit kompaktem Träger supp(Φ) ⊂ Ω.

Ist K ⊂ Ω kompakt, so definiert man

D_{K}(Ω) := {Φ ∈ X : supp(Φ) ⊂ K}

Mit der Seminorm

p_{K,m}(Φ) := sup_{|α|≤m_{m∈K}} |D^{α}Φ|

und der Nullumgebungsbasis

U(ϵ_{1}, …, ϵ_{n}; m_{1}, …, m_{n}; K) = {Φ ∈ D_{K}(Ω) : p_{K,m_{j}}(Φ) ≤ ϵ_{j} , 1 ≤ j ≤ n}

wird D_{K}(Ω) := X_{K} zum lokalkonvexen topologischen Vektorraum.

Definition 5.8 Mit der induktiven Topologie τ bzgl. {X_{K}} ist (X, τ) eine lokalkonvexer topologischer Raum, der mit D(Ω) bezeichnet wird. Der duale Raum D′(Ω) zu D(Ω) besteht aus allen stetigen linearen Funktionalen auf D(Ω). Ein Element T ∈ D′(Ω) heißt Distribution, eine Element ϕ ∈ D(Ω) heißt Testfunktion.

Definition 5.9 Sei T ∈ D′(Ω), α ein Multiindex. Die Ableitung D^{α}T ist die Distribution

(D^{α}T)(ϕ) = (−1)^{|α|}T (D^{α}ϕ) .

Sei nun u ∈ L_{1}^{loc}(Ω), d.h. u|_{A} ∈ L_{1}(A) für jede meßbare Teilmenge A ⊂ Ω, für die Ā ⊂ Ω und Ā kompakt ist. Sei außerdem

T_{u}ϕ := ∫_{Ω} u(x)ϕ(x)dx .

Das Funktional T_{u} ist nach Definition linear und sogar stetig, d.h. T_{u} ist Distribution bzw. T_{u} ∈ D′(Ω).

Definition 5.10 Sei u ∈ L_{1}^{loc}(Ω). Die Funktion v ∈ L_{1}^{loc}(Ω) heißt schwache Ableitung D^{α}u von u, falls

T_{v} = D^{α}(T_{u}) ,

d.h.

∫_{Ω} v(x)ϕ(x)dx = (−1)^{|α|} ∫ u(x)D^{α}ϕ(x)dx für alle ϕ ∈ D(Ω).

Dieses führt nun zur gesuchten Definition von Sobolev-Räumen:

Definition 5.11 (Sobolev-Räume über schwache Ableitungen) Sei m ∈ ℕ_{0} , p ∈ ℝ , p ≥ 1 , Ω ⊂ ℝ^{n}. W^{m,p}(Ω) sei die Menge aller Funktionen u, die schwache Ableitungen D^{α}u ∈ L^{p}(Ω) besitzen. Man setzt

∥u∥_{m,p} = [ ∑_{|α|≤m}(∥D^{α}u∥_{p})^{p}]^{1/p} ,

womit W^{m,p}(Ω) zum Banachraum wird.

Um sich das Verhältnis zwischen W^{m,p}(Ω) und H^{m,p}(Ω) verständlich zumachen, betrachtet man den Raum X := C^{∞}(Ω) ∩ H^{m,p}(Ω). X ist dicht in H^{m,p}(Ω). Jedes Element x ∈ X liegt auch in W^{m,p}(Ω). Es folgt:

H^{m,p}(Ω) ⊂ W^{m,p}(Ω)

Für 1 ≤ p < ∞ gilt sogar

H^{m,p}(Ω) = W^{m,p}(Ω).

Für die Konvergenzuntersuchungen bei Randwertaufgaben zweiter bzw. vierter Ordnung sind vor allem die Sobolev-Räume H^{1,2}(Ω) = W^{1,2}(Ω) und H^{2,2}(Ω) = W^{2,2}(Ω) von Bedeutung. Der Raum H^{m,2}(Ω) bildet mit dem Skalarprodukt

(u, v) := ∫_{Ω}( ∑_{|α|≤m} D^{α}uD^{α}v)dx

einen Hilbert-Raum, d.h. einen vollständigen normierten Raum mit Skalarprodukt, der alle im weiteren verlangten Eigenschaften bieten wird.

5.2 Variationsgleichungen

5.2.1 Beispiel 1: Lineares Randwertproblem

Als erstes und einfachstes Beispiel einer Variationsgleichung betrachten wir zunächst das lineare Randwertproblem

−u′′(x) + d(x)u′(x) + c(x)u(x) = f(x) in Ω := (0, 1)  (21)

mit u(0) = u(1) = 0 .

Dabei seien d, c und f gegebene stetige Funktionen und das Problem besitze eine Lösung, d.h. es existiere eine Funktion u ∈ C^{2}(Ω), die der Differentialgleichung genügt. Damit gilt für diese Lösung u trivialerweise auch

∫_{Ω}(−u′′ + du′ + cu)φ dx = ∫_{Ω} fφ dx

für beliebige Funktionen φ ∈ C_{0}^{∞}(Ω). Genügt umgekehrt eine Funktion u ∈ C^{2}(Ω) der Integralgleichung, so genügt u auch der zugrunde gelegten Differentialgleichung. Integriert man nun partiell, so erhält man zuerst

−u′φ|_{x=0}^{1} + ∫_{Ω} u′φ′ dx + ∫_{Ω}(du′ + cu)φ dx = ∫_{Ω} fφ dx

und unter Verwendung der Randwerte dann schließlich

∫_{Ω} u′φ′ dx + ∫_{Ω}(du′ + cu)φ dx = ∫_{Ω} fφ dx.

Also hat die zugehörige Variationsgleichung die Form

a(u, φ) = F (φ) für alle φ ∈ C_{0}^{∞}(Ω) ,  (22)

wobei

a(u, φ) := ∫_{Ω}[u′φ′ + (du′ + cu)φ] dx ,
F (φ) := ∫_{Ω} fφ dx .

5.2.2 Beispiel 2: Poissonsche Differentialgleichung

Als nächstes Beispiel einer einfachen elliptischen Randwertaufgabe bietet sich die Poissonsche Differentialgleichung mit homogenen Dirichlet-Randwerten an:

−Δ u(x_{1}, x_{2}) = f(x_{1}, x_{2}) in Ω ⊂ ℝ^{2}  (23)

mit u|_{Γ} = 0 , wobei Γ =: ∂Ω stückweise glatt ist.

Mit entsprechenden Funktionen φ ∈ C_{0}^{∞}(Ω) erhält man analog zum ersten Beispiel mittels partieller Integration eine zugehörige schwache Variante:

− ∫_{Ω} Δuφ dx = ∫_{Ω} fφ dx

⇒ −(\frac{∂u}{∂x_{1}} + \frac{∂u}{∂x_{2}} )φ|_{Γ} + ∫_{Ω} (\frac{∂u}{∂x_{1}} \frac{∂φ}{∂x_{1}} + \frac{∂u}{∂x_{2}} \frac{∂φ}{∂x_{2}} ) dx = ∫_{Ω} fφ dx

⇒ ∫_{Ω} (\frac{∂u}{∂x_{1}} \frac{∂φ}{∂x_{1}} + \frac{∂u}{∂x_{2}} \frac{∂φ}{∂x_{2}} ) dx = ∫_{Ω} fφ dx
⇒ a(u, φ) = F (φ) für alle φ ∈ C_{0}^{∞}(Ω) mit φ|_{Γ} = 0  (24)

wobei

a(u, φ) := ∫_{Ω} (\frac{∂u}{∂x_{1}} \frac{∂φ}{∂x_{1}} + \frac{∂u}{∂x_{2}} \frac{∂φ}{∂x_{2}} ) dx ,
F (φ) := ∫_{Ω} fφ dx .

5.2.3 Ausblick auf Lösungsmethoden

Erweitert man nun die Definitionsbereiche der Abbildungen a und F auf einen geigneten Hilbert-Raum H, so lassen sich mit der zugehörigen Norm Aussagen formulieren, die die Existenz und Eindeutigkeit einer Lösung u ∈ H mit

a(u, φ) = F (φ) für alle φ ∈ H

sichern. Zur Lösung eines Variationsproblems kann man das klassische Galerkin-Verfahren anwenden, dessen Grundidee sich folgendermaßen darstellt (vgl.: [10], S. 116):

1. Wähle H_{n} = span{h_{1}, …, h_{n}} ⊂ H.

2. Bestimme u_{n} ∈ H_{n} mit a(u_{n}, φ_{n}) = F (φ_{n}) für alle φ_{n} ∈ H_{n}.

Es kommt also darauf an, endlichdimensionale Untervektorräume H_{n} von H zu finden, die einerseits geometrisch einfach sind und andererseits gute Approximationseigenschaften besitzen. Diese Eigenschaften werden besonders gut von den Räumen finiter Elemente erfüllt. Hierbei wird das Gebiet Ω in möglichst gleichmäßige Polyeder, z.B. Dreiecke oder Quadrate im zweidimensionalen Fall, zerlegt (bzw. durch eine solche Zerlegung approximiert, falls der Rand von Ω gekrümmt ist). Finite Elemente sind dann stückweise polynomiale Funktionen von vorgegebenem Grad, im einfachsten Fall also stückweise lineare Funktionen auf Dreiecken. (vgl.: [20], S. 171)

Die weitere Vorgehensweise besteht dann darin, die Basisfunktionen h_{k} so zu wählen, daß jedes h_{k}

1. einen kleinen kompakten Träger besitzt,

2. stückweise glatt ist und

3. durch einfache Funktionen, z.B. Polynome, definiert ist.

Die Approximationsgüte hängt dann einerseits von der Polynomordnung und andererseits von der Feinheit der gewählten Gebietszerlegung ab. Meist ist es sinnvoll, die Polynomordnung festzuhalten, also z. B. nur stückweise lineare oder quadratische Elemente zuzulassen, und dafür die Zerlegung genügend zu verfeinern.

5.3 Das Variationsproblem im Hilbert-Raum

Die Existenz und Eindeutigkeit der Lösung u einer Variationsgleichung wird in einem geeigneten Hilbert-Raum, z. B. in dem Sobolev-Raum H^{m,p}(Ω), durch den Satz von Lax-Milgram gesichert, der einige Bedingungen an die Abbildungen a und F stellt. Diese sollen zuvor definiert werden:

Definition 5.12 Sei H ein Hilbert-Raum mit zugehörigem Skalarprodukt (·, ·) und zugehöriger Norm ∥·∥, a eine Bilinearform H×H → ℝ und F ein lineares Funktional H → ℝ. F heißt beschränkt, wenn es eine Konstante C gibt, so daß

∥F x∥ ≤ C∥x∥ für alle x ∈ H ,

und die Norm von F ist definiert durch

∥F ∥ := sup_{x∈H,x≠0} \frac{∥F x∥}{∥x∥} .

a heißt beschränkt, wenn es eine Konstante K gibt, so daß

a(u, v) ≤ K ∥u∥ ∥v∥ für alle u, v ∈ H ,

und die Norm von a ist definiert durch

∥a∥ := sup \frac{a(u, v)}{∥u∥ · ∥v∥} .

a heißt symmetrisch, falls für alle u, v ∈ H a(u, v) = a(v, u) gilt und a heißt koerziv, falls es eine Konstante α > 0 gibt, so daß

a(u, u) ≥ α∥u∥^{2} für alle u ∈ H .

Von zentraler Bedeutung für den Beweis den Satzes von Lax-Milgram ist der Rieszsche Darstellungssatz:

Satz 5.1 (Rieszscher Darstellungssatz) Jedes lineare beschränkte Funktional F auf einem Hilbert-Raum H ist von der Form

F (φ) = (φ, b) ,

wobei b ∈ H durch F eindeutig bestimmt ist. Umgekehrt definiert (φ, b) für jedes b ∈ H ein beschränktes, lineares Funktional auf H. Es gilt ∥F ∥ = ∥b∥.

(zum Beweis siehe [10], S. 120)

5.3.1 Satz von Lax-Milgram

Satz 5.2 (Lax-Milgram) Sei a eine beschränkte, koerzive Bilinearform auf dem Hilbert-Raum H und F ein beschränktes lineares Funktional auf H. Dann gibt es ein eindeutiges Element u ∈ H, so daß

a(φ, u) = F (φ) für alle φ ∈ H .

Ferner gilt die Abschätzung ∥u∥ ≤ \frac{1}{α}∥F ∥ für ein α > 0.

Beweis: Sei v ∈ H. Die Abbildung

φ ⟼ a(φ, v)

ist ein beschränktes lineares Funktional auf H. Aufgrund des Rieszschen Satzes existiert ein eindeutiges u_{v} ∈ H mit

a(φ, v) = (φ, u_{v}) für alle φ ∈ H .

Die Abbildung T : H → H mit v ↦ u_{v} ist beschränkt und linear:

• Beschränktheit: Nach Definition gilt

∥u_{v}∥^{2} = (u_{v}, u_{v}) = a(u_{v}, v) ≤ ∥a∥ · ∥u_{v}∥ · ∥v∥ ,

so daß mit der Beschränktheit von a folgt:

∥T v∥ = ∥u_{v}∥ ≤ ∥a∥ · ∥v∥ ≤ K · ∥v∥

für eine Konstante K ∈ ℝ.

• Linearität: Für v_{1}, v_{2} ∈ H gilt

a(φ, v_{1}) = (φ, T v_{1}) und a(φ, v_{2}) = (φ, T v_{2}) .

Dann gilt aufgrund der Bilinearität von a

a(φ, v_{1} + v_{2}) = (φ, T v_{1} + T v_{2}) ,
d.h. T (v_{1} + v_{2}) = T v_{1} + T v_{2} .

Durch die Koerzvität von a folgt schließlich noch, daß

α∥v∥^{2} ≤ a(v, v) = (v, T v) ≤ ∥v∥ · ∥T v∥ ,
⇒ ∥T v∥ ≥ α∥v∥ .

T ist also injektiv und T^{−1} ist beschränkt; der Wertebereich R(T) ist abgeschlossen. Damit folgt bereits R(T) = H, d.h. T ist surjektiv. Sonst gäbe es aufgrund des Projektionssatzes (siehe z.B. [10], S. 119) mit M = R(T) ein Element

z ∈ M^{⊥} ⊂ H , z ≠ 0 .

D.h. für alle v ∈ H würde a(z, v) = (z, T v) gelten, insbesondere für z = v:

a(z, z) = (z, T z) = 0 ⇒ z = 0 ,

was den erwünschten Widerspruch zur Voraussetzung liefert. Die inverse Abbildung T^{−1} ist deshalb eine lineare, beschränkte Abbildung von H auf H. Also gilt:

a(φ, T^{−1}w) = (φ, w) für alle φ, w ∈ H .

Die Anwendung des Rieszschen Satzes auf die Abbildung F liefert ein eindeutiges Element ũ ∈ H, für das gilt:

F (φ) = (φ, ũ) für alle φ ∈ H .

Es folgt für eben dieses ũ ∈ H und für alle φ ∈ H

F (φ) = (φ, ũ) = a(φ, T^{−1}ũ) .

Setzt man nun u := T^{−1}ũ, so hat man eine eindeutige Lösung u der Variationsgleichung erhalten. (q.e.d.)

5.3.2 Anwendung auf ein Modellproblem

Wir betrachten erneut die Poissongleichung mit homogenen Dirichlet Randbedingungen bzw. die zugehörige Variationsgleichung (siehe Beispiel 2): a(u, φ) = F (φ), wobei

a(u, φ) := ∫_{Ω} (\frac{∂u}{∂x_{1}} \frac{∂φ}{∂x_{1}} + \frac{∂u}{∂x_{2}} \frac{∂φ}{∂x_{2}} ) dx ,
F (φ) := ∫_{Ω} fφ dx ,

In diesem Beispiel ist der geeignete Hilbert-Raum H = H_{0}^{1,2}(Ω) mit der Norm

∥u∥_{1,2} := [∫_{Ω} ∑_{|α|≤1} |D^{α}u(t)|^{2} dt]^{1/2} ,

welche äquivalent ist zu der Halbnorm

|u|_{1,2} := [∫_{Ω} ∑_{|α|=1} |D^{α}u(t)|^{2} dt]^{1/2} ,

wie die Poincarésche Ungleichung zeigt:

Satz 5.3 (Poincarésche Ungleichung) Sei Ω ⊂ ℝ^{n} , d > 0. Ω liege zwischen den Ebenen x_{n} = 0 und x_{n} = d (d.h. Ω ist Lebesgue-meßbar). Dann gibt es eine Konstante K mit

|u|_{m,p} ≤ ∥u∥_{m,p} ≤ K · |u|_{m,p} für alle u ∈ H_{0}^{m,p}(Ω)

(Zum Beweis siehe z.B. [20], S. 159.)

Für das betrachtete Modellproblem der Poissongleichung liefert uns das

a(u, u) = ∫_{Ω} (\frac{∂u}{∂t_{1}} )^{2} + (\frac{∂u}{∂t_{2}} )^{2} dt
= ∫_{Ω} ∑_{|α|=1} |D^{α}u|^{2} dt
= (|u|_{1,2})^{2}

und damit die geforderte Koerzivität der Abbildung a:

a(u, u) = (|u|_{1,2})^{2} ≥ λ · (∥u∥_{1,2})^{2}

für ein λ > 0. Die nächste geforderte Eigenschaft erhält man aus der Hölderschen Ungleichung: Satz 5.4 (Höldersche Ungleichung) Ist u ∈ L^{p}(Ω) , v ∈ L^{q}(Ω) , 1/p + 1/q = 1, so gilt

∫_{Ω} u v ≤ ∥u∥_{p} · ∥v∥_{q}

(Zum Beweis siehe z.B. [20], S. 278.) Wendet man diese Ungleichung auf a an, erhält man die Beschränktheit von a:

a(u, φ) = ∫_{Ω} Du · Dφ dt
≤ ∥Du∥_{2} · ∥Dφ∥_{2}
= |u|_{1,2} · |φ|_{1,2}
⇒ a(u, φ) ≤ C · ∥u∥_{1,2} · ∥φ∥_{1,2} für eine Konstante C ∈ ℝ .

Da die Bilinearität von a offensichtlich ist, sind alle Anforderungen an a erfüllt und es bleibt nun noch die Abbildung F zu betrachten. Für f ∈ L^{2}(Ω) stellt

F : H_{0}^{1,2}(Ω) → ℝ mit φ ↦ ∫_{Ω} fφ dx

ein beschränktes lineares Funktional dar. Es ist nämlich

∥F ∥ := sup_{φ≠0} \frac{|F φ|}{∥φ∥} ≤ ∥f∥_{2} < ∞ ,

denn nach der Hölderschen und der Poincaréschen Ungleichung gilt:

F φ = ∫_{Ω} fφ ≤ ∥f∥_{2} · ∥φ∥_{2} ≤ ∥f∥_{2} · ∥φ∥_{1,2} .

Also sind für f ∈ L^{2}(Ω) alle Voraussetzungen für den Satz von Lax-Milgram gegeben und es existiert somit eine eindeutige Lösung des Variationsproblems der Poissonschen Differentialgleichung.

5.3.3 Das Optimierungsproblem

Die Variationsgleichung a(u, φ) = F (φ) kann auch als Optimierungsproblem gestellt werden, wie der folgende Satz zeigt:

Satz 5.5 Sei a eine beschränkte, symmetrische, koerzive Bilinearform auf dem Hilbert-Raum H, F (φ) := (b, φ) mit b ∈ H und J definiert als

J(φ) := \frac{1}{2}a(φ, φ) − F (φ) .

Es gilt: Ein Element u ∈ H löst das Optimierungsproblem

J(φ) ⟶ min ! bei φ ∈ H

genau dann, wenn u der Variationsgleichung a(u, φ) = F (φ) genügt. Es existiert genau ein u ∈ H, welches J(φ) minimiert.

Beweis: Sei zunächst u Lösung der Variationsgleichung und φ ∈ H. Nach Definition gilt

J(φ) − J(u) = \frac{1}{2}a(φ, φ) − F (φ) − \frac{1}{2}a(u, u) + F (u)
= \frac{1}{2}a(φ − u, φ − u) + a(u, φ − u) − F (φ − u)
= \frac{1}{2}a(φ − u, φ − u) ≥ 0
⇒ J(φ) ≥ J(u) für alle φ ∈ H ,  (25)

also löst u das Optimierungsproblem.

Genügt nun umgekehrt u dem Optimierungsproblem, so folgt

J(u) ≤ J(u + λφ) für alle λ ∈ ℝ , φ ∈ H

und damit nach Definition

\frac{1}{2}a(u, u) − F (u) − \frac{1}{2}a(u + λφ, u + λφ) + F (u + λφ) ≤ 0
⇒ \frac{1}{2}λ^{2}a(φ, φ) + λa(u, φ) − λF (φ) ≥ 0
⇒ \frac{1}{2}λa(φ, φ) + a(u, φ) − F (φ) ≥ 0 .  (26)

Betrachtet man nun den Grenzübergang von λ ↓ 0, so erhält man, daß einerseits

a(u, φ) − F (φ) ≥ 0 ,

und andererseits für λ ↑ 0, daß

a(u, φ) − F (φ) ≤ 0 .

Also ist u Lösung der Variationsgleichung a(u, φ) = F (φ). Die Eindeutigkeit dieser Lösung folgt direkt aus dem Satz von Lax-Milgram. (q.e.d.)

6 Diskretisierung der Probleme

Wie im vorherigen Abschnitt gezeigt, ist die Lösung einer klassischen Randwertaufgabe auch Lösung der schwachen Formulierung bzw. Variationsgleichung. Umgekehrt liefern Lösungen der Variationsgleichung unter zusätzlichen Regularitätseigenschaften (entsprechende Glattheit) auch Lösungen des passenden Randwertproblems. D. h. man löst die Problemstellung in geeigneten Hilbert-Räumen, wie H_{0}^{m,p}(Ω), um die Lösung einem solchen Raum zuordnen zu können und so Aussagen über die Glattheit der Lösung zu gewinnen.

Da nur in den seltensten Fällen eine Lösung explizit durch Integration o. ä. berechnet werden kann, ist man auf approximative numerische Verfahren angwiesen. Diese können nur in endlichen Dimensionen arbeiten, so daß man gezwungen ist, Probleme in unendlich dimensionalen Räumen durch solche in endlich dimensionalen zu beschreiben.

Die im Folgenden beschriebenen numerischen Techniken basieren auf der approximativen Lösung der Variationsgleichung

a(u, φ) = F (φ) für alle φ ∈ H  (27)

bzw. des Optimierungsproblems

J(φ) = \frac{1}{2}a(φ, φ) − F (φ) ⟶ min ! bei φ ∈ H ,  (28)

indem anstelle des zugrunde liegenden unendlich dimensionalen Raumes H ein Teilraum

H_{n} ⊂ H

mit dim H_{n} < ∞ gewählt wird. Statt des ursprünglichen Problems sind dann Aufgaben der Form

J(φ_{n}) = \frac{1}{2}a(φ_{n}, φ_{n}) − F (φ_{n}) ⟶ min ! bei φ_{n} ∈ H_{n}  (29)

zu lösen. Da wegen H_{n} ⊂ H auch H_{n} mit dem gleichen Skalarprodukt (·, ·) einen Hilbert-Raum bildet und a(·, ·) über H_{n} ebensolche Eigenschaften wie über H besitzt, läßt sich die abstrakte Theorie auch auf das endlich dimensionale Problem anwenden. Das diskretisierte Problem besitzt somit nach Lax-Milgram eine eindeutige Lösung u_{n} ∈ H_{n}. Diese läßt sich also auch durch die notwendige und hinreichende Bedingung

a(u_{n}, φ_{n}) = F (φ_{n}) für alle φ_{n} ∈ H_{n}  (30)

charakterisieren. Das auf dem diskreten Optimierungsproblem (29) basierende Verfahren heißt Ritz-Verfahren, das auf der diskreten Variationsgleichung (30) basierende Verfahren wird Galerkin-Verfahren genannt. Auf das Ritz-Verfahren wird in weiteren nicht näher eingegangen, da das Galerkin-Verfahren für das hier betrachtete Problem geeigneter ist.

6.1 Galerkin-Verfahren

Zur Konvergenz des Galerkin-Verfahrens gilt der häufig als Lemma von Cea bezeichnete folgende Satz:

Satz 6.1 (Cea) Sei a(·, ·) eine beschränkte, koerzive Bilinearform. Dann sind für jedes beschränkte, lineare Funktional F die Aufgaben (27) und (30) eindeutig lösbar. Für die zugehörigen Lösungen u ∈ H und u_{n} ∈ H_{n} gilt die Abschätzung

∥u − u_{n}∥ ≤ \frac{M}{γ} · inf_{φ_{n}∈H_{n}} ∥u − φ_{n}∥ ,

für zwei Konstanten M , γ ∈ ℝ.

Beweis: Die eindeutige Lösbarkeit von (30) folgt mit H_{n} anstelle von H aus dem Satz von Lax-Milgram, da sich wegen H_{n} ⊂ H die vorausgesetzten Eigenschaften von a auf H_{n} übertragen. Mit H_{n} ⊂ H folgt aus (27) auch

a(u, φ_{n}) = F (φ_{n}) für alle φ_{n} ∈ H_{n}

und unter Beachtung der Linearität erhält man mit (30) hieraus die Beziehung

a(u − u_{n}, φ_{n}) = 0 für alle φ_{n} ∈ H_{n} .

Insbesondere gilt also auch

a(u − u_{n}, u_{n}) = 0
⇒ a(u − u_{n}, u − u_{n}) = a(u − u_{n}, u − φ_{n}) für alle φ_{n} ∈ H_{n} .

Aus der Beschränktheit und Koerzivität von a folgt letztendlich

γ∥u−u_{n}∥^{2} ≤ M∥u−u_{n}∥·∥u−φ_{n}∥ für alle φ_{n} ∈ H_{n} und geeignete Konstanten γ , M ,

was die Behauptung liefert. (q.e.d.)

Da die Dimension von H_{n} endlich ist, gibt es endlich viele linear unabhängige h_{i}, i = 1, …, n, die den Teilraum H_{n} aufspannen, d.h.

H_{n} = span{h_{1}, …, h_{n}} = {v : v(x) = ∑_{i=1}^{n} s_{i}h_{i}(x)} ,

mit rellen Koeffizienten s_{i} für i = 1, …, n. Mit der Linearität von a(·, ·) und F ist damit (30) äquivalent zu

a(u_{n}, h_{i}) = F (h_{i}) , i = 1, …, n .

Beachtet man, daß die gesuchte Lösung u_{n} von (30) wegen u_{n} ∈ H_{n} eine Darstellung der Form

u_{n}(x) = ∑_{j=1}^{n} s_{j}h_{j}(x) , x ∈ Ω

besitzt, dann bildet dies ein lineares Gleichungssystem

∑_{j=1}^{n} a(h_{j}, h_{i})s_{j} = F (h_{i}) , i = 1, …, n  (31)

zur Bestimmung der Koeffizienten s_{j} ∈ ℝ , j = 1, …, n. Die Gleichungen (31) werden Galerkin-Gleichungen genannt. Eine wichtige Eigenschaft dieser Gleichungen läßt sich direkt aus den Eigenschaften von a ableiten:

Lemma 6.2 Ist a(·, ·) eine koerzive, beschränkte Bilinearform, so ist die Koeffizientenmatrix des Gleichungssystems (31) regulär.

Beweis: Sei z ∈ ℝ^{n} eine Lösung des homogenen Systems

∑_{j=1}^{n} a(h_{j}, h_{i})z_{j} = 0 , i = 1, …, n ,  (32)

dann gilt auch

∑_{i=1}^{n} ∑_{j=1}^{n} a(h_{j}, h_{i})z_{j}z_{i} = 0 .

Mit der Bilinearität von a erhält man hieraus

a(∑_{i=1}^{n} h_{i}z_{i}, ∑_{j=1}^{n} h_{j}z_{j}) = 0

und wegen der Koerzivität dann schließlich

∑_{i=1}^{n} h_{i}z_{i} = 0 .

Da die Funktionen h_{i} als linear unabhängig voraus gesetzt wurden, muß bereits z = 0 gelten. Also hat das homogene System (32) nur die triviale Lösung. Damit ist die Koeffizientenmatrix (im weiteren als Steifigkeitsmatrix bezeichnet) von (30) reglär. (q.e.d.)

Insgesamt ergibt sich die folgende schematische Übersicht zu den Variationsgleichungen und ihren diskreten Formulierungen:

Galerkin-Gleichungen ∑_{j=1}^{n} a(h_{j}, h_{i})s_{j} = F (h_{j})

i = 1, …, n

Galerkin-Verfahren a(u_{n}, φ_{n}) = F (φ_{n}) für alle φ_{n} ∈ H_{n} Ritz-Verfahren J(φ_{n}) ⟶ min! bei φ_{n} ∈ H_{n} Variationsgleichung a(u, φ) = F (φ)

für alle φ ∈ H

Optimierungsproblem J(φ) ⟶ min!

bei φ ∈ H

⟺

bei Symmetrie

⟺

bei Symmetrie

[figure]

6.2 Erweiterungen auf nichtlineare Probleme

Die im vorherigen Abschnitt mit Hilfe des Satzes von Lax-Milgram gemachten Aussagen zur Existenz von Lösungen abstrakter Variationsgleichungen setzen voraus, daß es sich bei der betrachteten Abbildung a(·, ·) um eine Bilinearform und bei F um ein lineares Funktional handelt. Die bisher gegebenen Resultate beziehen sich also nur auf lineare Randwertprobleme. Unter gewissen Bedingungen lassen sich diese aber auf nicht lineare Differentialoperatoren verallgemeinern, indem z. B. kontraktive Operatoren in Verbindung mit dem Banachschen Fixpunktsatz oder auch monotone Iterationsschemata mit zugehörigen Kompaktheitsargumenten verwendet werden.

Es bezeichne wieder H einen Hilbert-Raum mit dem zugehörigen Skalarprodukt (·, ·) und B : H ⟶ H sei ein Operator mit den folgenden Eigenschaften:

1. Es existiert ein γ > 0 derart, daß

(Bu − Bφ, u − φ) ≥ γ∥u − φ∥^{2}

für alle u, φ ∈ H.

2. Es existiert ein M > 0 derart, daß

∥Bu − Bφ∥ ≤ M∥u − φ∥

für alle u, φ ∈ H.

Als abstraktes Ausgangsproblem wird die Operatorgleichung

Bu = 0  (33)

untersucht. Diese ist äquivalent zur nichtlinearen Variationsgleichung

(Bu , φ) = 0 für alle φ ∈ H .  (34)

In direkter Verallgemeinerung des Satzes von Lax-Milgram gilt

Satz 6.3 Unter den Voraussetzungen 1. und 2. besitzt die Operatorgleichung (33) eine eindeutige Lösung u ∈ H. Diese ist Fixpunkt des durch

T_{r}φ := φ − rBφ , φ ∈ H ,

definierten Operators T_{r} : H ⟶ H, der für Parameterwerte r ∈ (0, \frac{2γ}{M^{2}} ) kontrahierend ist.

Beweis: Wir zeigen zunächst die Kontraktivität des eingeführten Operators T_{r} für die angebebenen Parameterwerte. Nach Konstruktion gilt mit der durch das Skalarprodukt erzeugten Norm ∥u∥ := \sqrt{(u, u)} des Hilbert-Raumes

∥T_{r}y − T_{r}φ∥^{2} = ∥y − rBy − [φ − rBφ]∥^{2}
= (y − φ − r[By − Bφ] , y − φ − r[By − Bφ])
= ∥y − φ∥^{2} − 2r(By − Bφ, y − φ) + r^{2}∥By − Bφ∥^{2}
≤ (1 − 2γr + r^{2}M^{2})∥y − φ∥^{2}

für alle y, φ ∈ H. Also ist T_{r} kontrahierend für

|r^{2}M^{2} − 2γr + 1| < 1
⟺ r ∈ (0, \frac{2γ}{M^{2}} ) .

T_{r} besitzt damit nach dem Banachschen Fixpunktsatz einen eindeutigen Fixpunkt u ∈ H, d. h. es gilt

u = T_{r}u = u − rBu .

Folglich löst u auch die Operatorgleichung (33), so daß nun nur noch die Eindeutigkeit der Lösung zu zeigen bleibt. Es sei also ũ ∈ H ebenfalls Lösung der Operatorgleichung, d. h.

Bũ = Bu = 0 .

Mit der starken Monotonieeigenschaft (1., siehe oben) des Operators und der Bilinearität des Skalarproduktes folgt, daß

0 = (Bu − Bũ, u − ũ) ≥ γ∥u − ũ∥^{2} ≥ 0 .

Somit gilt ũ = u, so daß die Lösung von (33) eindeutig ist. (q.e.d.)

Ebenso lassen sich Gleichungen mit Operatoren A : H^{∗} ⟶ H anstelle von B betrachten. Dabei bezeichnet H^{∗} den zu H gehörigen Dualraum, d.h. den Raum der beschränkten, linearen Funktionale F : H ⟶ ℝ. Es sei vorausgesetzt, daß A den folgenden Eigenschaften genüge:

1. Der Operator A ist stark monoton, d. h. es existiert ein γ > 0 derart, daß

(Au − Aφ, u − φ) ≥ γ∥u − φ∥^{2} für alle u, φ ∈ H .

2. Der Operator A ist Lipschitz-stetig, d. h. es existiert eine Konstante M > 0

derart, daß

∥Au − Aφ∥ ≤ M∥u − φ∥ für alle u, φ ∈ H .

Dann besitzt das Problem

Au = F  (35)

für beliebiges F ∈ H^{∗} eine eindeutige Lösung u ∈ H. Diese Aussage erhält man unmittelbar aus Satz 6.3 mit Hilfe des für φ ∈ H durch

Bφ := T (Aφ − F)

definierten Operators B : H ⟶ H. Analog zum Beweis des Satzes von Lax-Milgram bezeichnet T : H^{∗} ⟶ H den geeigneten Operator zum Rieszschen Darstellungssatz, der jedem beschränkten linearen Funktional F ∈ H^{∗} eine Element T (F) ∈ H mit

F (φ) = (T F, φ) für alle φ ∈ H

zuordnet. Das Problem (35) ist äquivalent zu

(Au, φ) = F (φ)  (36)

für alle φ ∈ H.

Unter den getroffenen Voraussetzungen kann nun auch das Lemma von Cea auf die betrachteten Aufgaben übertragen werden. Hierzu gilt:

Satz 6.4 Es seien A : H^{∗} ⟶ h stark monoton bezüglich einer Konstanten γ > 0 und Lipschitz-stetig bezüglich einer Konstanten M > 0 sowie F ∈ H^{∗}. Zu jedem linearen Unterraum H_{n} ⊂ H mit dim H_{n} < ∞ existiert ein eindeutig bestimmtes u_{n} ∈ H_{n}, das der diskreten Variationsgleichung

(Au_{n}, φ_{n}) = F (φ_{n}) für alle φ_{n} ∈ H_{n}  (37)

genügt. Mit der exakten Lösung u ∈ H der Variationsgleichung gilt hier zusätzlich die Abschätzung

∥u − u_{n}∥ ≤ \frac{M}{γ} · inf_{φ_{n}∈H_{n}} ∥u − φ_{n}∥ .

(Zum Beweis siehe [17], S. 112.)

Im Unterschied zu linearen Randwertaufgaben erhält man im vorliegenden Fall dann natürlich auch nichtlineare Galerkin-Gleichungen. Wird der Ansatz

u_{n}(x) = ∑_{j=1}^{n} s_{j}h_{j}(x)

gewählt, so ist (37) äquivalent zum Gleichungssystem

(A(∑_{j=1}^{n} s_{j}h_{j}) , h_{i}) = F (h_{i}) , i = 1, …, n .

Auf die iterativen Methoden, die zur Lösung solcher Systeme nötig sind, wird später noch eingegangen.

6.3 Methode der finiten Elemente

6.3.1 Zusammenfassung

In diesem Abschnitt soll in Anlehnung an [32], [26] und [17] die zugrundeliegende Idee und das sich daraus ergebende Vorgehen der Methode der finiten Elemente zur Lösung von Aufgaben, wie sie in den vorherigen Abschnitten eingeführt worden sind, beschrieben werden.

Wie in Abschnitt 6.2.3 bereits skizziert, besitzt die Methode der finiten Elemente die folgenden vier typischen Merkmale:

1. Zerlegung des Grungebietes Ω ⊂ ℝ^{n} in geometrisch einfache Teilgebiete Ω_{n};

2. Definition von geeigneten Ansatzfunktionen u_{j} und linear unabhängigen, stückweise stetigen Formfunktionen N_{j} über den Teilgebieten;

3. iterative Lösung des erhaltenen Gleichungssystemes liefert Approximationen u_{n};

4. Auswertung des Fehlers zwischen exakter Lösung u und u_{n} bezüglich geeigneter Norm.

Zu 1. Die gegebene Aufgabe wird diskretisiert, indem das Grundgebiet Ω in einfache Teilgebiete, die sogenannten Elemente Ω_{n}, zerlegt wird. Im Fall von zweidimensionalen Problemen wird das Grundgebiet in Dreiecke, Parallelogramme, krummlinige Dreiecke oder Vierecke eingeteilt. Selbst wenn nur geradlinige Elemente verwendet werden, erreicht man mit einer entsprechend feinen Diskretisierung eine gute Approximation des Grundgebietes, wobei krummlinige Elemente selbstverständlich die Güte der Annäherung erhöhen. Bei räumlichen Problemen erfolgt eine Diskretisierung des dreidimensionalen Gebietes in Tetraederelemente, Quaderelemente oder anderen, dem Problem angepaßten, möglicherweise auch krummflächig berandeten Elementen. Die wichtigsten Finite-Elemente-Ansätze im ℝ^{2} und ℝ^{3} findet man z.B. in [17], S. 129ff, abgebildet und ausführlich beschrieben.

Zu 2. In jedem der Elemente wird für die gesuchte Funktion ein problemgerechter Ansatz bzw. problemgerechte Formfunktionen N_{j} gewählt, um so eine Approximation der gesuchten Funktion u bzw. eine Ansatzfunktion u_{n} für u der Form

u_{n} = ∑_{j=1}^{n} s_{j} · N_{j}

mit unbekannten Koeffizienten s_{j} zu finden.

Die Art des Ansatzes hängt dabei einerseits von der Form des Elementes ab; andererseits kann auch das zu behandelnde Problem den zu wählenden Ansatz beeinflussen, da die Ansatzfunktionen beim Übergang von einem Element in ein benachbartes ganz bestimmte problemabhängige Stetigkeits- und Differenzierbarkeitsbedingungen und Randwerte erfüllen müssen. Elemente mit Ansatzfunktionen, welche den Stetigkeitsbedingungen genügen, heißen konform.

Um diese Anforderungen tatsächlich zu erfüllen, ist der Funktionsverlauf innerhalb eines betrachteten Elementes durch gegebene Funktionswerte oder auch durch Werte partieller Ableitungen in bestimmten Punkten des Elementes, den Knotenpunkten, auszudrücken. Mit Hilfe solcher Werte als Koeffizienten stellt sich die Ansatzfunktion als Linearkombination der Formfunktionen dar. Werden ausschließlich vorgegebene Funktionswerte zur Bestimmung der Ansatzfunktionen verwendet, so bezeichnet man die zugehörigen Elemente als Lagrange-Elemente. Im Unterschied dazu nutzen Hermite-Elemente dem Prinzip der Hermiteschen Interpolation folgend auch Vorgaben von Ableitungswerten zur Bestimmung der Ansatzfunktionen.

Nimmt man konkret den Fall an, daß in den m Knotenpunkten p_{j}^{(e)} des Elementes e nur Funktionswerte u(p_{j}^{(e)}) verwendet werden, so erhält die Ansatzfunktion die Darstellung

u_{n}^{(e)} = ∑_{j=1}^{m} u(p_{j}^{(e)}) · N_{j}^{(e)} .

Da diese Darstellung für den Funktionswert u(p_{j}^{(e)}) jedes Knotenpunktes gelten muß, so muß die Formfunktion N_{j}^{(e)} notwendigerweise im Punkt p_{j}^{(e)} gleich 1 sein und in den anderen Knotenpunkten des Elementes (e) verschwinden. Um an dieser Stelle die Grundlage für die Anwendung des Galerkin-Verfahrens im Sinne der Methode der finiten Elemente vorzubereiten, betrachtet man die globale Darstellung der gesuchten Funktion u_{n} im ganzen Grundgebiet, bestehend aus der Vereinigungsmenge der Elemente. Numeriert man die Werte in den Knotenpunkten fortlaufend von 1 bis n, dann läßt sich das Ergebnis der Zusammensetzung formulieren als

u_{n} = ∑_{j=1}^{n} u(p_{j}) · N_{j} .

Daraus wird ersichtlich, daß die globalen Formfunktionen N_{j} nur in denjenigen Elementen von Null verschieden sind, welche den Knotenpunkt p_{j} gemeinsam haben, so daß also die Funktionen N_{j} nur in einem sehr beschränkten Teilgebiet von Null verschieden sind, also nur einen lokalen Träger aufweisen, was ein wesentliches Charakteristikum der Methode der finiten Elemente ist. Die Formfunktionen müssen also zwei Eigenschaften erfüllen:

• N_{i} hat polynomiale Gestalt über jedem Element und ist eindeutig festgelegt durch die Funktionswerte an den Knotenpunkten dieses Elementes.

• N_{i}(p_{j}) = δ_{ij}, wobei δ_{ij} das Kronecker-Symbol darstellt, das für i = j 1 und sonst 0 ist.

In einem Ansatz dieser Form lassen sich Randbedingungen einfach durch Vorgabe entsprechender Werte für die betreffenden Knotenpunkte verwirklichen.

Zu 3. und 4. soll in den nächsten Abschnitten genauer auf einige iterative Methoden und auf die Methode der gewichteten Residuen eingegangen werden.

6.3.2 Methode der gewichteten Residuen

Das Ziel der oben beschriebenen Methoden ist es also, eine approximative Lösung u_{n} der Form

u_{n} = ∑_{j=1}^{n} s_{j} · N_{j}

mit unbekannten Koeffizienten s_{j} zu finden. Wird nun dieser Ansatz in eine partielle Differentialgleichung der Form

L (u) = 0

mit einem geeignet definierten Differentialoperator L eingesetzt, so wird sie in den seltensten Fällen erfüllt sein, vielmehr resultiert ein sogenanntes Residuum R = L (u_{n}). Da der genaue Fehler ebenso wie die exakte Lösung meist unbekannt sind, verlangt man also, daß das Residuum im Inneren des Gebietes Ω möglichst klein werden soll. Die Idee, die hinter der Methode der Gewichteten Residuen steht, ist nun folgende:

Gesucht sind Koeffizienten s_{1}, …, s_{n} derart, daß das gewichtete Mittel des Residuums R über Ω verschwindet für n linear unabhängige Gewichtsfunktionen oder auch Testfunktionen W_{i}:

∫_{Ω} R · W_{i} = 0 , i = 1, …, n

In der Methode von Galerkin werden die n Gewichtsfunktionen der Reihe nach gleich den n Formfunktionen N_{i} gewählt, welche die Bedingung der linearen Unabhängigkeit ja bereits erfüllen. Mit dieser Wahl der Gewichtsfunktionen erreicht man, daß die Residuenfunktion orthogonal zum Funktionenunterraum ist, der von den N_{1}, …, N_{n} aufgespannt wird. Diese Tatsache rechtfertigt das beschriebene Vorgehen, indem die resultierende approximative Lösung u_{n} in diesem Sinne die bestmögliche im Raum der Ansatzfunktionen darstellt.

6.3.3 Modellproblem: Die 1D Poisson-Gleichung

Man betrachte das eindimensionale Randwertproblem

−u′′(x) = f(x) , für x ∈ Ω = [0, 1] , u(0) = u_{L} , u(1) = u_{R} .

Das Gebiet Ω = [0, 1] wird unterteilt in m Elemente der Länge h.

Ω_{e} = [(e − 1)h , eh] mit h = 1/m und e = 1, …, m .

Die Knotenpunkte seien bezeichnet mit p_{i} = (i − 1)h , i = 1, …, n = m + 1.

Jede Formfunktion N_{i} soll nun stückweise linear sein mit N_{i}(p_{j}) = δ_{ij}, was es recht einfach macht, die vorgegebenen Randwerte zu erfüllen. Wähle dazu die Funktion

ψ(x) = u_{L}N_{1}(x) + u_{r}N_{n}(x) ,

die folgenden Ansatz liefert:

u(x) ≈ u_{n}(x) = ψ(x) + ∑_{j=2}^{n−1} s_{j}N_{j}(x) = ∑_{j=1}^{n} s_{j}N_{j}(x) ,

wobei natürlich s_{1} = u_{L} und s_{n} = u_{R} sein müssen.

Die Gleichung der gewichteten Residuen erhält man nun, wie schon beschrieben, durch Einsetzen der Darstellung von u_{n} in die Differentialgleichung, Multiplikation mit den Gewichtsfunktionen W_{i} = N_{i} , i = 1, …, n und anschließender Integration über das ganze Intervall [0, 1]. Vergleichbar mit dem Variationsproblem der Poissonschen Differentialgleichung in Abschnitt 5.2.2 erhält man mittels partieller Integration folgende diskrete Formulierung:

∑_{j=1}^{n}( ∫_{0}^{1} N′_{i}N′_{j}) s_{j} = ∫_{0}^{1} fN_{i} für i = 1, …, n .

Dieses Gleichungssystem läßt sich mit folgenden Bezeichnungen in Matrixschreibweise darstellen:

A := (a_{ij})_{i,j=1}^{n} mit a_{ij} := ∫_{0}^{1} N′_{i}N′_{j} ,
b := (b_{i})_{i=1}^{n} mit b_{i} := ∫_{0}^{1} fN_{i} ,
s := (s_{i})_{i=1}^{n} .

Man erhält also das Gleichungssystem

A s = b .

Dabei wird A wie schon zuvor erwähnt in Anlehnung an Probleme der Elastitätstheorie Steifigkeitsmatrix genannt.

Eine sehr feine Zerlegung des Gebietes Ω führt zwar dazu, daß die erzeugte Matrix A eine große Dimension besitzt, allerdings ist diese Matrix dann nur schwach besetzt, da die vorher aufgezeigten Eigenschaften der Formfunktionen N_{i} zu diversen Nulleinträgen führen: Der Matrix-Eintrag a_{ij} ist genau dann von Null verschieden, wenn die Knotenpunkte p_{i} und p_{j} zum selben Element gehören. In anderen Worten:

a_{ij} ≠ 0 ⟺ j ∈ {i − 1, i, i + 1} .

Setzt man die betrachteten Knotenpunkte auf den Rand eines jeden Elementes fest, so könnten die Formfunktionen wie in Abbildung 13 konstruiert sein. Nimmt man nun der Einfachheit halber eine konstante Elementlänge bzw. -größe h an, so folgt direkt, daß N′_{i} = ± 1/h. Betrachtet man nun die rechte Seite von Abbildung 13, so erhält man für die Randeinträge und weiteren Nicht-Null-Einträge von A:

a_{1,1} = \frac{1}{h} ,
a_{i,i−1} = ∫_{0}^{1} N′_{i−1}N′_{i} = −\frac{1}{h} für i = 2, …, n ,
a_{i,i} = ∫_{0}^{1} N′_{i}N′_{i} = \frac{2}{h} für i = 2, …, n − 1 ,
a_{i,i+1} = ∫_{0}^{1} N′_{i+1}N′_{i} = −\frac{1}{h} für i = 1, …, n − 1 ,
a_{n,n} = \frac{1}{h} .

Zur genauen Betrachtung der Ausdrücke b_{i} verwendet man die explizite Darstellung der N_{i} mit

N_{i}(x) = \{ \begin{matrix}0  , & x ≤ p_{i−1} \\ (x − p_{i−1})/h , & p_{i−1} < x ≤ p_{i} \\ (p_{i+1} − x)/h , & p_{i} < x ≤ p_{i+1} \\ 0  , & x > p_{i+1} \end{matrix}

für 1 < i < n. Die Formeln für N_{1} und N_{n} müssen leicht abgeändert werden. Die Integrale der Ausdrücke b_{i} können nur in den seltensten Fällen analytisch berechnet werden, so daß häufig ein numerisches Integrationsschema nach der Trapezregel angewendet wird. Hierbei werden die Knotenpunkte auch zu den Integrationspunkten:

b_{i} = ∫_{0}^{1} f N_{i} ≈ \frac{h}{2}f(p_{1})N_{i}(p_{1}) + ∑_{j=2}^{n−1} h f(p_{j})N_{i}(p_{j}) + \frac{h}{2}f(p_{n})N_{i}(p_{n}) .

1  2  3  4  5  6  7  x  1  2  3  4  5  6  7  x

Abbildung 13: Skizze der stückweise linearen Formfunktionen N_{i} in einer Dimension für einige Elemente. Die Skalierung an der x-Achse bezieht sich auf die Knotenpunkte des Gitters. Die Graphen links illustrieren also die Formfunktionen N_{3}, N_{4}, N_{5} und N_{6}, während die Graphen rechts die jeweiligen Ableitungen skizzieren. (aus [26], S. 126 und 128)

Das liefert b_{1} = f(0)·h/2, b_{i} = f(p_{i})·h für i = 2, …, n−1 und schließlich auch b_{n} = f(1)·h/2. Die Randwerte kann man allerdings auch direkt mit s_{1} = u_{L} und s_{n} = u_{R} in das Gleichungssystem einfließen lassen. Also hat das globale Gleichungssystem die folgende Gestalt, die formale Ähnlichkeit mit der finiten Differenzen-Methode aufweist:

A s = b

⟺ \{ \begin{matrix}s_{1} & = & u_{L} \\ −\frac{1}{h}s_{i−1} + \frac{2}{h}s_{i} − \frac{1}{h}s_{i+1} & = & f(p_{i})h für i = 2, …, n − 1 \\ s_{n} & = & u_{r}x \end{matrix}

mit p_{i} = (i − 1)h , i = 2, …, n.

6.3.4 Elementmatrizen

Im vorherigen Abschnitt wurde das zu lösende Gleichungssystem direkt hergeleitet, was für mehrdimensionale Probleme mit geometrisch anspruchsvolleren Gebieten nicht durchführbar ist.

Um in diesen Fällen das globale lineare Gleichungssystem aufzustellen, betrachtet man zunächst nur die jeweiligen Elementmatrizen und setzt diese zur globalen Matrix zusammen. Man schreibt

a_{ij} = ∫_{0}^{1} N′_{i}N′_{j} dx = ∑_{e=1}^{m} a_{ij}^{(e)} , wobei a_{ij}^{(e)} = ∫_{Ω_{e}} N′_{i}N′_{j}dx ,

und b_{i} = ∫_{0}^{1} fN_{i} dx = ∑_{e=1}^{m} b_{i}^{(e)} , wobei b_{i}^{(e)} = ∫_{Ω_{e}} fN_{i}dx, .

Nun ist, wie schon zuvor, a_{ij}^{(e)} genau dann ungleich Null, wenn p_{i} und p_{j} Knotenpunkte des Elementes Ω_{e} sind, da sonst die Formfunktionen N′_{i} N′_{j} verschwinden. Genau so ist auch b_{i}^{(e)} nur dann ungleich Null, wenn p_{i} Knotenpunkt von Ω_{e} ist. Es ist also sinnvoll, nur die lokalen Knotenpunkte des gerade behandelten Elementes zur Berechnung von a_{ij}^{(e)} und b_{i}^{(e)} heran zu ziehen, was zur Notwendigkeit einer lokalen Knotennummerierung für jedes Element führt. In unserem einfachen linearen Fall mit den Elementen Ω_{e} = [p_{e}, p_{e+1}] erhält p_{e} die lokale Knotennummer 1 und dementsprechend p_{e+1} die lokale Knotennummer 2.

Es gibt also eine eindeutige Zuordnung i = q(e, r), die zu jeder lokalen Knotennummer r des Elementes e die zugehörige globale Knotennummer i liefert und umgekehrt. Im betrachteten linearen Fall ist

q ( e , r) = e − 1 + r für r = 1, 2 .

Bei mehrdimensionalen Problemen läßt sich diese Zuordnung nicht mehr so einfach analytisch angeben, sondern existiert nur in Form einer Zuordnungstabelle für alle Knotenpunkte. Man sammelt also im betrachteten Fall zunächst die Nicht-Null-Einträge von A_{ij}^{(e)} in einer 2 × 2 Elementmatrix Ã_{rs}^{(e)}, wobei r, s = 1, 2 die lokalen Knotennummern bezeichnen, und ebenso die Elementvektoren b̃_{r}^{(e)} für r = 1, 2. Vor allem bei mehrdimensionalen Problem ist es dann noch zwingend notwendig, jedes pysikalische Element vor weiteren Berechnungen auf ein Referenz-Element mit fester Einheitsgröße zu transformieren. Dieses kann mit einer vollständigen Variablentransformation bei den Knotenpunkten, den Formfunktionen und sämtlichen Integralen und Ableitungen erreicht werden. Siehe dazu z.B. [26], S. 132-136.

Abbildung 14: Zusammensetzung der globalen Matrix: Die zu Element Nr. 3 und Nr. 4 gehörigen Elementmatrizen werden zur globalen Matrix hinzu gefügt. (nach [26], S. 137)

Zum Schluß bleibt dann nur noch die Aufgabe, die einzelnen Elementmatrizen und Elementvektoren zu einem globale Gleichungssystem zusammenzusetzen, was mittels der lokalen Knotennummern (r, s) und den Zuordnungen q(e, r) = i, q(e, s) = j auf einfach Weise realisiert werden kann. Abbildung 14 zeigt, wie die 2 × 2 Elementmatrizen unseres linearen Modellproblemes zu einer globalen Matrix zusammengefügt werden. Vollendet man diesen Prozess, so erhält man das erwartete globale Gleichungssystem für konstante Elementlänge h:

( \begin{matrix}\frac{1}{h} & −\frac{1}{h} & 0 & ⋯ & ⋯ & ⋯ & 0 \\ −\frac{1}{h} & \frac{2}{h} & −\frac{1}{h} & 0 & ... & ... & ⋮ \\ 0 & −\frac{1}{h} & \frac{2}{h} & −\frac{1}{h} & ... & ... & ⋮ \\ ⋮ & ... & ... & ... & ... & ... & ⋮ \\ ⋮ & ... & ... & ... & ... & ... & 0 \\ ⋮ & ... & ... & ... & −\frac{1}{h} & \frac{2}{h} & −\frac{1}{h} \\ 0 & ⋯ & ⋯ & ⋯ & 0 & −\frac{1}{h} & \frac{1}{h} \end{matrix} ) ( s_{1} ; s_{2} ; ⋮ ; ⋮ ; ⋮ ; s_{n−1} ; s_{n} ) = ( u_{L} ; f(p_{2})h ; ⋮ ; ⋮ ; ⋮ ; f(p_{n−1})h ; u_{R} )

6.4 Iterative Methoden

Nachdem man partielle Differentialgleichungen mit der Finite-Elemente Methode diskretisiert hat, bleibt nun noch die Aufgabe, das wie oben erhaltene lineare Gleichungssystem effektiv und doch genau zu lösen. In diesem Abschnitt sei das Gleichungssystem ganz allgemein mit

A · x = b

für eine Matrix A ∈ ℝ^{(n×n)} und einen gegebenen Vektor b ∈ ℝ^{n} bezeichnet.

Lineare Systeme kann man einerseits direkt analytisch und andererseits iterativ lösen, wobei die direkte Lösung häufig das Eliminationsverfahren von Gauss beinhaltet. Im allgemeinen benötigt das Gaussverfahren im Laufe der Elimination also n^{2} Einträge der Matrix A und führt ungefähr 2/3 · n^{3} Rechenoperationen durch, d.h. Additionen, Subtraktionen, Multiplikationen und Divisionen. Gerade für mehrdimensionale Probleme und die daraus resultierenden sehr großen Gleichungssysteme sind die iterativen Methoden gegenüber den direkten Methoden klar im Vorteil: Iterative Methoden können aufgrund der schwachen Besetztheit vieler Koeffizienten-Matrizen deutlich effektiver arbeiten, da sie auch deutlich weniger Rechenoprationen benötigen. Meist basieren die iterativen Löser auf Matrix-Vektor-Multiplikationen, bei denen sich die vielen Null-Einträge dann positv auf den Rechenaufwand auswirken.

In diesem Zusammenhang ist das Maß für die Effektivität eines iterativen Lösers natürlich die Konvergenzgeschwindigkeit des Iterationsschemas. Um diese Konvergenzgeschwindigkeit zu erhöhen, gibt es verschiedenste Möglichkeiten der Vorkonditionierung eines Gleichungssystemes. Are Magnus Bruaset verdeutlicht in [5] besonders anschaulich den gravierenden Unterschied zwischen der Effektivität vorkonditionierter iterativer Methoden und direkter Methoden. Er betrachtet die Finite-Differenzen-Diskretisierung der 3D Poisson-Gleichung und vergleicht die Zeiten, die theoretisch zur Lösung des Gleichungssystemes benötigt werden. Vergleicht man die Werte der Gauss-Elimination mit denen der vorkonditionierten Methode der konjugierten Gradienten (s. [5], Tabelle 1.1, S. 5), so stellt man fest, das die Unterschiede mit zunehmender Anzahl der Unbekannten immer deutlicher werden. Bei 8 000 000 Unbekannten würde die iterative Methode theoretisch ca. 2 600 Sekunden benötigen, während die direkte Methode unpraktikable 832 Jahre beschäftigt wäre.

6.4.1 Klassische Methoden: Jacobi und Gauss-Seidel Verfahren

Im folgenden sollen in Anlehnung an die in der später noch genauer erläuterten C++-Biliothek Diffpack verwendeten Methoden (s. [26], S. 592ff) einige klassische Iterationsmethoden in ihren wesentlichen Zügen erläutert werden. Zu genaueren Details siehe [16], S. 147ff, und [5], S. 10ff.

Unter den klassischen Iterationmethoden versteht man meist Verfahren wie Jacobi, Gauss-Seidel und das SOR-Verfahren. Diese Verfahren basieren alle auf einer Zerlegung der Koeffizientenmatrix des Gleichungssystemes

A · x = b

für eine Matrix A ∈ ℝ^{(n×n)} und einen gegebenen Vektor b ∈ ℝ^{n}. Schreibt man A in der Form A = M − N, wobei M eine reguläre (n × n)-Matrix ist, so verändert sich das zugrunde liegende Gleichungssystem zu

M x = N x + b .

Diese Zerlegung der Koeffizientenmatrix wird in bezug auf die Matrix M mit ’matrix splitting’ oder Vorkonditionierung betitelt, wobei M dann als Vorkonditionierer bezeichnet wird. Die Regularität der Matrix M liefert, daß Systeme der Form Mv = w durch nur O(n) Rechenoperationen zu lösen sind. Also ist es auch keine Schwierigkeit, das Gleichungssystem durch Multiplikation mit M^{−1} noch weiter zu modifizieren:

x = M^{−1}N x + M^{−1}b := G x + c

Diese umformulierte Variante des Gleichungssystemes legt nun offensichtlich das folgende Iterationsschema nahe: Für gegebenen Anfangsvektor x^{0} berechne

x^{k} = G x^{k−1} + c für k = 1, 2, …

Die oben definierte Matrix G spielt bei den Konvergenzbetrachtungen der iterativen Methode eine wesentliche Rolle, wie man direkt an dem Zusammenhang

x^{k} − x = G^{k}(x^{0} − x)

erkennen kann. Eine iterative Methode heißt konvergent, falls die Folge {x^{k}} für k → ∞ und einen Startvektor x^{0} gegen die exakte Lösung x konvergiert. Man sieht also, daß die Bedingung

lim_{k→∞} ∥G^{k}∥ = 0

notwendig und hinreichend für die Konvergenz gegen x ist. Diese Bedingung ist äquivalent zur Formulierung

ρ(G) < 1

wobei mit ρ(G) der Spektralradius von G, also das Maximum der Beträge der Eigenwerte, bezeichnet ist. Zum Beweis dieser Äquivalenz siehe z.B. [17], S. 257, Lemma 5.2. Verkleinert man ρ(G), so erhöht sich die Konvergenzgeschwindigkeit.

Von besonderem Interesse ist dann außerdem die Konvergenzrate bzw. relative Konvergenzgeschwindigkeit R_{k}, die definiert ist durch

R_{k}(G) = −\frac{ln ∥G^{k}∥}{k} ,

und damit die zugehörige asymptotische Konvergenzrate liefert:

R_{∞}(G) = lim_{k→∞} R_{k}(G) = − ln ρ(G)

Um den Anfangsfehler zum Beispiel um einen Faktor ϵ zu verringern, d.h. um

∥x − x^{k−1}∥ ≤ ϵ · ∥x − x^{0}∥

zu erreichen, benötigt man − ln ϵ/R_{∞}(G) Iterationschritte.

Da die Schreibweise

x^{k} = G x^{k−1} + c

für die spätere Implementierung nicht geeignet ist, wählt man meist die etwas abgewandelte Schreibweise

x^{k} = x^{k−1} + M^{−1}r^{k−1} ,

wobei r^{k−1} = b − Ax^{k−1} das Residuum nach dem Iterationsschritt k −1 bezeichnet. Zerlegt man nun noch A in Diagonalmatrix D, obere Dreiecksmatrix U und untere Dreiecksmatrix L, also A = L+D+U, dann lassen sich im folgenden die Unterschiede zwischen den klassischen Verfahren mittels dieser Notation darstellen. Desweiteren seien mit A_{r,s} die Blockmatrizen der Matrix A bezeichnet, also

A = ( \begin{matrix}A_{1,1} & ⋯ & A_{1,ν} \\ ⋮ & & ⋮ \\ A_{ν,1} & ⋯ & A_{ν,ν} \end{matrix} ) , x = ( x_{1} ; ⋮ ; x_{ν} ) und b = ( b_{1} ; ⋮ ; b_{ν} ) ,

wobei A_{r,s} ∈ ℝ^{(n_{r}×n_{s})} und ∑_{r=1}^{ν} n_{r} = n sein muß.

Jacobi-Verfahren Wählt man M als die Diagonale von A, also M = D und N = −L − U, so stellt sich das Iterationsschema fogendermaßen dar:

x^{k} = x^{k−1} + D^{−1}r^{k−1} ,

Schreibt man dieses Schema explizit aus, so folgt für den k-ten Iterationschritt und die i-te Komponente der gesuchten Lösung x^{k}

x_{i}^{k} = x_{i}^{k−1} + \frac{1}{a_{i,i}} (b_{i} − ∑_{j=1}^{n} a_{i,j} · x_{j}^{k−1})

für alle i = 1, …, n. Zur Betrachtung des Konvergenzverhaltens wählt man wiederum die Darstellung x^{k} = G_{J} x^{k−1} + c, wobei bei diesem Verfahren

G_{J} = D^{−1}(L + U) und c = D^{−1}b

sind. Das führt nun zur Konvergenzbedingung ρ(D^{−1}(L + U)) < 1. Handelt es sich bei der Matrix A nun z.B. um eine positiv definite symmetrische Matrix, oder allgemeiner um eine Matrix mit Diagonaldominanz, so sieht man direkt, daß

ρ(G_{J}) = ρ(D^{−1}(L + U)) ≤ ∥D^{−1}(L + U)∥_{∞} = max_{i} ∑_{j≠i} |a_{i,j}/a_{i,i}| < 1

ist. Je ”stärker“ die Diagonaldominanz, desto höher die Konvergenzgeschwindigkeit.

Gauss-Seidel-Verfahren Ein lineares Gleichungssystem, dessen Koeffizientenmatrix eine obere oder untere Dreiecksgestalt besitzt, ist bekannterweise besonders einfach zu lösen, so daß der Gedanke, M als obere Dreicksmatrix zu wählen, nahe liegt. Es ist also M = D + L und demnach dann N = −U. Die Gauss-Seidel-Methode hat also allgemein folgende Gestalt:

(D + L) · x^{k} = −U x^{k−1} + b

Der Standard-Algorithmus zum Lösen eines Gleichungssystemes mit Koeffizientenmatrix in oberer Dreicksgestalt angewandt auf diese Formel liefert für den k-ten Iterationschritt und die i-te Komponente der gesuchten Lösung x^{k}

x_{i}^{k} = \frac{1}{a_{i,i}} (b_{i} − ∑_{j=1}^{i−1} a_{i,j}x_{j}^{k} − ∑_{j=i+1}^{n} a_{i,j}x_{j}^{k−1})

für alle i = 1, …, n. Um Aussagen bez. des Konvergenzverhaltens machen zu können, betrachtet man wiederum die Darstellung x^{k} = G_{GS} x^{k−1} + c, wobei hier

G_{GS} = (D − L)^{−1} · U und c = (D − L)^{−1} · b

sind. Das führt zur Konvergenzbedingung ρ(G_{GS}) = ρ((D − L)^{−1} · U) < 1. Man kann nun für eine positv definite Koeffizientenmatrix A direkt Konvergenz des Verfahrens zeigen (siehe dazu z.B. [17], S. 259) oder noch allgemeiner die Gleichungungen

ρ(G_{GS}) = (ρ(G_{J}))^{2} und R_{∞}(G_{GS}) = 2 · R_{∞}(G_{J})

für jede Koeffizientenmatrix A mit Nicht-Null Einträgen auf der Diagonalen. Siehe dazu z.B. [5], S. 14.

Der wesentliche Unterschied der expliziten Iterationsschemata dieser beiden Verfahren liegt also in der Verwendung bereits berechneter Werte: Das Gauss-Seidel-Verfahren benutzt die neu berechnete Näherung bereits beim nächsten Teilschritt und wird aufgrund dessen auch als Einzelschritt-Verfahren bezeichnet, während das Jacobi-Verfahren die neu berechnete Näherung erst nach Vollendung aller Teilschritte einbringt und somit die Bezeichnung Gesamtschritt-Verfahren erhält.

6.4.2 Relaxationsverfahren: SOR und SSOR

Ein erheblicher Nachteil einfacher Iterationsverfahren vom Typ der Gauss-Seidel-und Jacobi-Verfahren besteht in der asymptotisch langsamen Konvergenz bzw. in der niedrigen asymptotischen Konvergenzrate R_{∞}. Gewisse Beschleunigungen dieser Verfahren lassen sich durch die zusätzliche Verwendung von Ralaxations-Parametern erreichen. Im weiteren soll ausschließlich das Gauss-Seidel-Verfahren als Grundmethode betrachtet werden.

SOR-Verfahren Bezeichnet x̃ einen durch einen Gauss-Seidel-Schritt vorhergesagten Wert für x im k-ten Iterationsschritt, so wird x^{k} als gewichtete Kombination aus x̃ und dem vorherigen Wert x^{k−1} dargestellt:

x^{k} = ωx̃ + (1 − ω)x^{k−1} .

Dabei ist ω > 0 der bereits oben genannte Relaxations-Parameter. Das Ziel ist, diesen Parameter so zu wählen, daß der Spektralradius ρ(G) der Iterationsmatrix minimal ist. Die resultierende Methode heißt ’successive overrelaxation’, kurz SOR, und stimmt für ω = 1 mit der Gauss-Seidel-Methode überein.

Die Relaxations-Methode kann ausgedrückt werden mittels

M = \frac{1}{ω} D + L und N = \frac{1 − ω}{ω} D − U .

Der zugehörige Algorithmus hat also offensichtlich die Form

x_{i}^{k} = \frac{ω}{a_{i,i}} (b_{i} − ∑_{j=1}^{i−1} a_{i,j}x_{j}^{k} − ∑_{j=i+1}^{n} a_{i,j}x_{j}^{k−1}) + (1 − ω)x_{i}^{k−1}

für i = 1, …, n.

SSOR-Verfahren Das SOR-Verfahren läßt sich durch Kombination von Vorwärtsund Rückwärtsvarianten symmetrisieren. Auf diese Weise erhält man das symmetrische SOR-Verfahren, kurz SSOR. Die Matrix M hat die Gestalt

M = \frac{1}{2 − ω} (\frac{1}{ω} D + L)(\frac{1}{ω} D)^{−1}(\frac{1}{ω} D + U) ,

so daß sich der zugehörige Algorithmus analog zum SOR-Verfahren erstellen läßt.

A. M. Bruaset vergleicht in [5], S. 18, die Konvergenzraten der vier erwähnten Iterationsmethoden in Anwendung auf ein Modellproblem: Betrachtet man die 2D Poisson-Gleichung −∇^{2}u = f auf einem mit Elementlänge h diskretisierten Einheitsrechteck mit u = 0 auf dem Rand, so erhält man folgende asymptotische Abschätzungen:

Konvergenzrate Jacobi Gauss-Seidel SOR SSOR

R_{∞}(G)  π^{2}h^{2}/2  π^{2}h^{2}  2πh > πh

Zu detaillierteren Erläuterungen der verwendeten Iterationsmethoden siehe [26], S. 591ff, und die Werke von A. Greenbaum [16] und A. M. Bruaset [5] zu iterativen Methoden und deren Analysis.

7 Existenz- und Eindeutigkeitssätze zu den Navier-Stokes Gleichungen

Viele Aussagen aus der Theorie zur Lösung partieller Differentialgleichungen lassen sich nicht auf die vollständigen Navier-Stokes-Gleichungen übertragen, da es sich hierbei um einen äußerst komplizierten instationären, nicht linearen Spezialfall handelt. Um trotzdem Existenz- und Eindeutigkeitssätze formulieren zu können, betrachtet man schwache bzw. Variationsformulierungen des Problems und versucht, mittels Konvergenz approximativer Verfahren und Abschätzungen geeigneter Hilbertraumnormen Existenz und Eindeutigkeit dieser schwachen Lösungen nachzuweisen. Im folgenden werden die klassischen Ergebnisse von Roger Temam [35] aus dem Jahre 1977 mit den aktuellen Ergebnissen von Prof. M. Wiegner [39] verglichen bzw. ergänzt, um ein vollständiges Bild des derzeitigen Wissensstandes zu erhalten.

7.1 Klassische Resultate: R. Temam, 1977

Das Problem der vollen Navier-Stokes-Gleichungen ist das folgende:

Gesucht ist eine vektorielle Funktion u : Ω × [0, T] → ℝ^{n} und eine skalare Funktion p : Ω × [0, T] → ℝ derart, daß

\frac{∂ u}{∂ t} + (u · ∇) u = −grad p + \frac{1}{Re} Δ u + f in Ω_{t} = Ω × (0, T)
div u = 0 in Ω_{t}
u = 0 auf ∂Ω × (0, T)
u(x, 0) = u_{0}(x) in Ω .

(38)

7.1.1 Notationen

f und u_{0} sind gegeben und auf ganz Ω ⊂ ℝ^{n} bzw. Ω_{t} ⊂ ℝ^{n+1} definiert. Zur Vereinfachung sei im weiteren \frac{1}{Re} = ν = 1 angenommen.

Die bis auf weiteres betrachteten Räume sind

H = {u ∈ C_{0}^{∞}(Ω) : div u = 0} ,
H = Vervollständigung von H in H_{0}^{1,2}(Ω) ,
L = Vervollständigung von H in L^{2}(Ω) .

Die Dualräume von H und L, d.h. die Räume der linearen Funktionale auf H bzw. L seien mit H^{∗} und L^{∗} bezeichnet. (·, ·) sei das zum Lebesgue-Raum L gehörige Skalarprodukt

(u, v) = ∫_{Ω} u · v ,

((·, ·)) das zum Hilbertraum H gehörige Skalarprodukt

((u, v)) = ∫_{Ω} [∑_{i=1}^{n}(D_{i}u · D_{i}v)] .

Zu den verwendeten Normen vergleiche Kapitel 5 und [35].

7.1.2 Schwache Lösungen des Linearen Problems

Das lineare Navier-Stokes-Problem hat folgende Gestalt:

\frac{∂ u}{∂ t} − Δ u + grad p = f in Ω_{t}
div u = 0 in Ω_{t}
u = 0 auf ∂Ω × [0, T]
u(x, 0) = u_{0}(x) in Ω

(39)

mit gegebenen Funktionen f und u_{0}.

Sind nun u und p klassische Lösungen von (39), also u ∈ C^{2}(Ω) , p ∈ C^{1}(Ω), so liefert wie schon zuvor das Produkt mit einer beliebigen Funktion φ ∈ H in Ω und anschließende Integration über Ω

∫_{Ω} \frac{∂ u}{∂ t} · φ − ∫_{Ω} Δu · φ + ∫_{Ω} grad p · φ = ∫_{Ω} fφ .  (40)

Die einzelnen Summanden kann man auf folgende Weise umformen:

∫_{Ω} \frac{∂ u}{∂ t} · φ = (\frac{∂ u}{∂ t}, φ) = (u_{t}, φ)

Da dies für alle φ ∈ H gilt, kann man dieses Skalarprodukt auch schreiben als

(u_{t}, φ) = \frac{d}{dt}(u, φ) .

Analog liefert einerseits

− ∫_{Ω} Δu · φ = ∇u · φ|_{∂Ω} + ∫_{Ω} ∇u∇φ

= 0 + ∫_{Ω} [∑_{i=1}^{n}(D_{i}u · D_{i}v)]

= ((u, φ))

und andererseits

∫_{Ω} fφ = (f, φ) ,

während

∫_{Ω} grad p · φ = p · φ|_{∂Ω} − ∫_{Ω} p · div φ = 0

Eingesetzt in (40) erhält man also die folgende schwache Formulierung des Problems (39):

Finde für gegebenes f ∈ L^{2}(0, T ; H^{∗}) und u_{0} ∈ L ein u ∈ L^{2}(0, T ; H), welches der Variationsgleichung

\frac{d}{dt}(u, φ) + ((u, φ)) = (f, φ) .  (41)
mit u(0) = u_{0}

für alle φ ∈ H genügt.

Satz 7.1 Für gegebenes f ∈ L^{2}(0, T ; H^{∗}) und u_{0} ∈ L existiert eine eindeutige Funktion u ∈ L^{2}(0, T ; H), die das Problem (41) erfüllt.

Beweis Man benutzt zum Beweis der Existenz die Galerkin-Methode. H ist separabel, d.h es existiert in H eine abzählbare dichte Teilmenge von Teilräumen H_{n} = span{h_{1}, …, h_{n}} mit linear unabhängigen Elementen h_{i} ∈ H. Also existiert eine in H vollständige Folge (h_{n})_{n∈ℕ} mit linear unabhängigen Folgegliedern.

Für jedes n ∈ ℕ definiert man nun eine approximative Lösung u_{n} von (41) wie folgt:

u_{n} = ∑_{i=1}^{n} s_{in}(t) · h_{i}  (42)

und mit \frac{∂ u}{∂ t} = u′

(u′_{n}, h_{j}) + ((u_{n}, h_{j})) = (f, h_{j}) , j = 1, …, n ,
u_{n}(0) = u_{0n} ,  (43)

wobei u_{0n} ∈ H_{n} mit u_{0n} ⟶ u_{0} für n → ∞.

Die Funktionen s_{in}(t) , 1 ≤ i ≤ n, sind skalare Funktionen definiert auf [0, T], zu deren Bestimmumg man mittels (42) und (43) ein System von Differentialgleichungen erhält.

∑_{i=1}^{n}(h_{i}, h_{j})s′_{in}(t) + ∑_{i=1}^{n}((h_{i}, h_{j}))s_{in}(t) = (f, h_{j}) , j = 1, …, , n;  (44)

Die Elemente h_{1}, …, h_{n} sind linear unabhängig, so daß direkt folgt, daß die Matrix mit den Einträgen (h_{i}, h_{j}) , 1 ≤ i, j ≤ n, regulär ist. Indem man nun diese Matrix invertiert, reduziert sich das Differentialgleichungssystem (44) auf ein lineares Gleichungssystem mit konstanten Koeffizienten.

s′_{in}(t) + ∑_{j=1}^{n} α_{ij}s_{jn}(t) = ∑_{j=1}^{n} β_{ij}(f(t), h_{j}) , i = 1, …, n ,  (45)

wobei α_{ij}, β_{ij} ∈ ℝ. Die Anfangsbedingung von (43) ist äquivalent zu den Gleichungen

s_{in}(0) = i-te Komponente von u_{0n} .  (46)

Das lineare Gleichungssystem (45) zusammen mit den Anfangsbedingungen (46) liefert nun eindeutige Lösungen für die s_{in} auf dem ganzen Intervall [0, T].

Da die skalaren Funktionen t ↦ (f, h_{j})(t) zweimal differenzierbar sind, sind es die Funktionen s_{in} auch, so daß für alle n ∈ ℕ gilt:

u_{n} ∈ L^{2}(0, T ; H) und u′_{n} ∈ L^{2}(0, T ; H) .  (47)

Als nächstes sollen unabhängig von n Abschätzungen für die Funktionen u_{n} gefunden werden, die eine Betrachtung des Grenzüberganges für n → ∞ erlauben. Dafür kehrt man wieder zurück zu den Ausgangsgleichungen

(u′_{n}, h_{j}) + ((u_{n}, h_{j})) = (f, h_{j}) , j = 1, …, n ,

multipliziert mit s_{jn} und summiert die Gleichungen für j = 1, …, n auf. Man erhält:

(u′_{n}(t), u_{n}(t)) + ∥u_{n}(t)∥^{2} = (f(t), u_{n}(t)) .

Aufgrund von (47) gilt:

2(u′_{n}(t), u_{n}(t)) = \frac{d}{dt}|u_{n}(t)|^{2}
⇒ \frac{d}{dt}|u_{n}(t)|^{2} + 2∥u_{n}(t)∥^{2} = 2(f(t), u_{n}(t))

≤ 2 · ∥f(t)∥ · ∥u_{n}(t)∥
≤ ∥u_{n}(t)∥^{2} + ∥f(t)∥^{2}
⇒ \frac{d}{dt}|u_{n}(t)|^{2} + ∥u_{n}(t)∥^{2} ≤ ∥f(t)∥^{2}  (48)
⇒ \frac{d}{dt}|u_{n}(t)|^{2} ≤ ∥f(t)∥^{2} .  (49)

Integriert man (49) von 0 bis s, T ≥ s > 0, so erhält man die Ungleichung:

|u_{n}(s)|^{2} ≤ |u_{0n}|^{2} + ∫_{0}^{s} ∥f(t)∥^{2} dt
≤ |u_{0}|^{2} + ∫_{0}^{T} ∥f(t)∥^{2} dt
⇒ sup_{s∈[0,T]} |u_{n}(s)|^{2} ≤ |u_{0}|^{2} + ∫_{0}^{T} ∥f(t)∥^{2} dt < ∞  (50)

Die rechte Seite von (50) ist also unabhängig von n endlich, d.h.:

Die Folge u_{n} verbleibt in einer beschränkten Menge von L^{∞}(0, T ; L) . (51)

Ebenso erhält man durch Integration von (48) von 0 bis T die Abschätzung

|u_{n}(T)|^{2} + ∫_{0}^{T} ∥u_{n}(t)∥^{2} dt ≤ |u_{0}|^{2} + ∫_{0}^{T} ∥f(t)∥^{2} dt ,

die damit zeigt:

Die Folge u_{n} verbleibt in einer beschränkten Menge von L^{2}(0, T ; H) . (52)

Die Bedingungen (51) und (52) ermöglichen nun die Betrachtung des Grenzüberganges für n → ∞:

Dank (51) existiert ein Element u ∈ L^{∞}(0, T ; L) und eine Teilfolge u_{n′} von u_{n}, derart, daß

u_{n′} ⟶ u unter der schwachen Topologie von L^{∞}(0, T ; L)

d.h. für jedes v ∈ L^{1}(0, T ; L) gilt:

∫_{0}^{T} (u_{n′}(t) − u(t) , v(t)) dt ⟶ 0 für n′ → ∞ .  (53)

Durch (52) erhält man zusätzlich, daß die Teilfolge u_{n′} in einem beschränkten Gebiet des L^{2}(0, T ; H) verbleibt. Das heißt ein erneuter Übergang zu einer wiedrum mit u_{n′} bezeichneten Teilfolge liefert die Existenz eines u_{∗} ∈ L^{2}(0, T ; H) derart, daß

u_{n′} ⟶ u_{∗} unter der schwachen Topologie von L^{2}(0, T ; H) ,

d. h., es gilt für alle v ∈ L^{2}(0, T ; H′):

∫_{0}^{T} (u_{n′}(t) − u_{∗}(t) , v(t)) dt ⟶ 0 für n′ → ∞ .

Vergleicht man dieses mit (52), so sieht man direkt, daß

∫_{0}^{T} (u(t) − u_{∗}(t) , v(t)) = 0

für alle v ∈ L^{2}(0, T ; L), also

⟹ u = u_{∗} ∈ L^{2}(0, T ; H) ∩ L^{∞}(0, T ; L) .

Um jetzt den Grenzübergang von (43) für n → ∞ zu betrachten, nimmt man zunächst skalare Funktionen ψ zur Hilfe, die stetig differenzierbar auf [0, T] sein sollen mit

ψ(T) = 0 .

Multipliziert man (43) mit ψ und integriert partiell, so erhält man

− ∫_{0}^{T} (u_{n} , h_{j} · ψ′) + ∫_{0}^{T} ((u_{n} , h_{j} · ψ)) = (u_{0n} , h_{j}) · ψ(0) + ∫_{0}^{T} (f , h_{j}) · ψ .

Betrachtet man nun an dieser Stelle den Grenzübergang für n = n′ ⟶ ∞, so gilt mit den zuvor gezeigten Abschätzungen:

− ∫_{0}^{T} (u , h_{j} · ψ′) + ∫_{0}^{T} ((u , h_{j} · ψ)) = (u_{0} , h_{j}) · ψ(0) + ∫_{0}^{T} (f , h_{j}) · ψ .

Summation über alle j = 1, …, n ergibt

− ∫_{0}^{T} (u , φ)ψ′ + ∫_{0}^{T} ((u , φ))ψ = (u_{0} , φ)ψ(0) + ∫_{0}^{T} (f , φ)ψ

für jedes φ, welches endliche Linearkombination der h_{j} ist.

Wählt man speziell ψ ∈ C_{0}^{∞}((0, T)), so folgt

\frac{d}{dt}(u, φ) + ((u, φ)) = (f, φ) für alle φ ∈ H .

Also ist die als Grenznzwert gewonnene Funktion u schwache Lösung des linearen Navier-Stokes Problemes. Zum Abschluß des Existenzbeweises bleibt nur noch zu zeigen, daß für dieses u auch u_{0} = u(0) gilt. Dazu multipliziert wiederum die schwache Fomulierung für diese Lösung u mit oben beschriebenen ψ, integriert partiell und erhält

− ∫_{0}^{T} (u, φ)ψ′ + ∫_{0}^{T} ((u, φ))ψ = (u(0), φ)ψ(0) + ∫_{0}^{T} (f, φ)ψ .

Vergleicht man dieses mit der aus der Grenzbetrachtung gewonnenen Gleichung, so folgt direkt, daß

(u(0) − u_{0} , φ) = 0 für alle φ ∈ H
⟹ u(0) = u_{0} ,

womit die Existenz einer Lösung bewiesen wäre. Zum Beweis der Stetigkeit und Eindeutigkeit siehe [35], S. 260-264.

7.2 Aktuelle Resultate: M.Wiegner, 1999

7.2.1 Der Stokes-Operator und die Stokes-Halbgruppe

Die analytische Behandlung partieller Differentialgleichungen basiert meist auf der Betrachtung spezifischer Differentialoperatoren in geeigneten Funktionenräumen. Bezieht man sich wie im letzten Abschnitt auf

H = {u ∈ C_{0}^{∞}(Ω) : div u = 0} ,
H = Vervollständigung von H in H_{0}^{1,2}(Ω) ,
L = Vervollständigung von H in L^{2}(Ω) ,

so stellt sich im Zusammenhang der Existenz- und Eindeutigkeitssätze die Frage: Gibt es eine eindeutige orthogonale Zerlegung

L^{2}(Ω) = L ⊕ G^{2}

mit einer linearen stetigen Abbildung

P : L^{2}(Ω) ⟶ L

und wie kann man G^{2} charakterisieren? Offensichtlich ist P ein von Ω abhängiger Operator, der bei Anwendung auf eine Differentialgleichung mit Lösung u weder die Randwerte von u erhält, noch mit Differentialoperatoren kommutiert, was in der Vergangenheit immer wieder zu falschen Behauptungen und Beweisen bez. der Navier-Stokes-Gleichungen geführt hat. Andererseits erhält P im allgemeinen die Stetigkeit und ist unabhängig vom Druck p. Es wird klar, daß das orthogonale Komplement zu den L-Funktionen mit div = 0 die Gradienten von L^{2}-Funktionen sein müssen. Wendet man den Operator P nun auf die Navier-Stokes-Gleichungen für f = 0 an, so erhält man

P (u_{t} − Δu + ∇p) = P (−u∇u) ,
wobei P u = u , da div u = 0 ,
P u_{t} = u_{t} ,
P (−u∇u) = −P (u∇u) ,
P (∇p) = 0 ,
so daß u_{t} \underbrace{−P (Δu)}_{=:Au} = −P (u∇u)
⇒ u_{t} + A u = −P (u∇u) .  (54)

Den durch Au = −P Δu definierten Operator bezeichnet man als Stokes-Operator oder Projektor. Mit der Hilfe dieses Operators kann das Ausgangs-Problem der partiellen Differentialgleichung auf das einfachere Problem einer gewöhnlichen Differentialgleichung zurückgeführt werden. Betrachtet man eine gewöhnliche Differentialgleichung der Form

y′ + A y = f

mit einer Matrix A ∈ Mat(n, ℝ) und Anfangswert y_{0}, so liefert die Methode der Variation der Konstanten (siehe [37], S. 25) eine Lösung der Form

y = e^{−tA}y_{0} − ∫_{0}^{t} e^{−(t−s)A}f ds .

Die Frage, ob man diese Eigenschaften von gewöhnlichen auf partielle Differentialgleichungen übertragen kann, und wie man e^{−tA} für einen Differentialoperator A verstehen kann, läßt sich durch die Halbgruppeneigenschaft des Stokes-Operators A erklären, d.h., daß

A infinitesimaler Erzeuger einer stetigen Halbgruppe e^{−tA}

ist. Diese Stokes-Halbgruppe ist beschränkt, was man z.B. über die Resolventenabschätzung gewinnen kann:

Satz 7.2 (Hille-Yosida) A : D(A) → L sei ein linearer Operator mit in dem Banachraum L ⊂ L^{2} dichtem Definitionsbereich D(A). Die Resolvente R(λ, A) = (λ I − A)^{−1} existiere für alle λ ∈ ℝ und es gelte

∥(I − \frac{1}{λ}A)^{−1}∥ ≤ 1 .

Dann erzeugt A eine eindeutig bestimmte, kontrahierende Halbgruppe.

(zum Beweis siehe [20], S. 133, und [13], S. 135)

Man hat also folgende Eigenschaften des Ausdruckes e^{−tA}:

1. Ist u_{0} ∈ L^{2}, so ist e^{−tA}u_{0} definiert als Lösung des Problems

u_{t} + A u = 0 , u(0) = u_{0} .

2. e^{−tA} : L → L ist beschränkt.

3. e^{−tA} erfüllt die Halbgruppeneigenschaft, d.h.

e^{−(t+s)A}u_{0} = e^{−tA}(e^{−sA}u_{0}) .

Es liegt also nahe, die Lösung der Navier-Stokes-Gleichungen durch eine Integralgleichung der Form

u = e^{−tA}u_{0} − ∫_{0}^{t} e^{−(t−s)A}P (u∇u) ds

zu konstruieren. Zum Beweis der Existenz der Lösung einer solchen Integralgleichung benötigt man zwei Abschätzungen der Stokes-Halbgruppe und die Höldersche Ungleichung für den Projektor P:

∥e^{−tA}v∥_{q} ≤ c_{0}t^{− \frac{n}{2} ( \frac{1}{p} − \frac{1}{q})}∥v∥_{p} für \frac{n}{2} ≤ p ≤ q < ∞ ,  (55)
∥∇e^{−tA}v∥_{n} ≤ c_{1}t^{− \frac{n}{2p}} ∥v∥_{p} für \frac{n}{2} < p ≤ n ,  (56)
∥P (v∇v)∥\frac{n}{1+δ} ≤ c_{2}∥v∥\frac{n}{δ} ∥∇v∥_{n} .  (57)

Vergleiche dazu [39], S. 5 und S. 9, und [22], S. 474.

7.2.2 Schwache Lösungen des Nichtlinearen Problems

Wie in Abschnitt 8.1 bereits gesehen, beschränken sich die klassischen Existenz- und Eindeutigkeitssätze auf schwache Lösungen der linearen Navier-Stokes Gleichungen. An dieser Stelle soll nun auch für die schwachen Lösungen des nichtlinearen Problems ein Existenz- und Eindeutigkeitsbeweis erbracht werden; siehe [39], S. 5-8. Die Notationen der betrachteten Funktionenräume werden aus dem vorherigen Abschnitt übernommen.

Eine schwache Lösung der Navier-Stokes-Gleichungen sei in folgendem Sinne definiert:

Definition 7.1 Sei f ∈ L^{2}(Ω_{t}) und u_{0} ∈ L. u ∈ L^{2}(0, T ; H) ist schwache Lösung der Navier-Stokes-Gleichungen, falls

∫_{0}^{T} ∫_{Ω}(−u · φ_{t} + ∇u∇φ + (u∇)u · φ) = ∫_{Ω} u_{0} · φ(0) + ∫_{0}^{T} ∫_{Ω} f · φ (58)

für alle φ ∈ H mit φ(T) = 0. Diese schwache Formulierung des Problems erhält man vergleichbar mit [35], S. 252, durch Multiplikation mit φ und anschließender Intergration über Ω_{t}.

Um die Eindeutigkeit einer solchen Lösung nachzuprüfen, bietet es sich an, die Gleichungen für zwei unterschiedliche Lösungen u_{1}, u_{2} zu subtrahieren, mit der Differenz v = u_{1} − u_{2} zu multiplizieren und anschließend über Ω_{t} = Ω × [0, T] zu integrieren. Diese Schritte liefern nun folgende Gleichung:

∫_{Ω_{t}} [(u_{1} − u_{2})_{t}(u_{1} − u_{2}) − (Δ(u_{1} − u_{2}))(u_{1} − u_{2}) + (u_{1}∇u_{1} − u_{2}∇u_{2})(u_{1} − u_{2})] = 0

(59)

Betrachtet man die einzelnen Integrale, so erhält man:

∫_{Ω_{t}}(u_{1} − u_{2})_{t}(u_{1} − u_{2}) = \frac{1}{2} ∫_{Ω} ∫_{0}^{T} (v^{2})_{t}
= \frac{1}{2} ∫_{Ω}[v^{2}(T) − v^{2}(0)]
= \frac{1}{2} ∫_{Ω} v^{2}(T) ,

− ∫_{Ω_{t}}(Δ(u_{1} − u_{2}))(u_{1} − u_{2}) = − ∫_{Ω_{t}} Δv · v
= − ∫_{0}^{T} ( \underbrace{∇v · v|_{∂Ω}}_{=0} − ∫_{Ω} ∇v · ∇v )

= ∫_{Ω_{t}}(∇v)^{2} ,

∫_{Ω_{t}} [(u_{1}∇)u_{1} − (u_{2}∇)u_{2}] (u_{1} − u_{2}) = ∫_{Ω_{t}} [(u_{1}∇)v + (v∇)u_{2}] v
= ∫_{Ω_{t}}(v∇)u_{2} · v + ∫_{Ω_{t}}(u_{1}∇)v · v
= ∫_{Ω_{t}}(v∇)u_{2} · v + \frac{1}{2} ∫_{Ω_{t}} u_{1}∇(v^{2})
= ∫_{Ω_{t}}(v∇)u_{2} · v + \frac{1}{2} ∫_{0}^{T} ( \underbrace{u_{1}v^{2}|_{∂Ω}}_{=0} − \underbrace{∫_{Ω} divu_{1} · v^{2}}_{=0} )

= ∫_{Ω_{t}}(v∇)u_{2} · v .

Also bleibt von Gleichung (59) nur noch folgendes zu betrachten :

\frac{1}{2} ∫_{Ω} v^{2}(T) + ∫_{Ω_{t}}(∇v^{2}) + ∫_{Ω_{t}}(v∇)u_{2} · v = 0

⟺ \frac{1}{2}∥v(T)∥_{2}^{2} + ∫_{0}^{T} ∥∇v∥_{2}^{2} = − ∫_{Ω_{t}}(v∇)u_{2} · v ≥ 0  (60)

Die rechte Seite schätzt man nun mittels Sobolev-Ungleichung und Youngscher Ungleichung ab. Die dafür notwendige Sobolev-Ungleichung ist

∥v∥_{4} ≤ c∥v∥_{2}^{1− \frac{n}{4}} ∥∇v∥2^{\frac{n}{4}} für n = 2, 3

und die Youngsche Ungleichung liefert zuerst

a · b ≤ ε · b^{p} + c_{ε} · a^{q}

für \frac{1}{p} + \frac{1}{q} = 1, und speziell für p = \frac{4}{n} , q = \frac{4}{4−n}

a · c^{\frac{4−n}{2}} · b^{\frac{n}{2}} ≤ ε · b^{2} + c_{ε} · a^{\frac{4}{4−n}} · c^{2} ,

so daß die folgenden Abschätzungen gerechtfertigt sind.

|∫_{Ω_{t}}(v∇)u_{2} · v| = |∫_{0}^{T} ∫_{Ω}((u_{1} − u_{2})∇)u_{2} · (u_{1} − u_{2})|
≤ ∫_{0}^{T} (∫_{Ω} |∇u_{2}|^{2})^{\frac{1}{2}} · ( ∫_{Ω} |(u_{1} − u_{2})^{2}|^{2})^{\frac{1}{2}}

= ∫_{0}^{T} ∥∇u_{2}∥_{2} · ∥u_{1} − u_{2}∥_{4}^{2}
(Sobolev Ungl.) ≤ ∫_{0}^{T} ∥∇u_{2}∥_{2} · (c∥u_{1} − u_{2}∥_{2}^{2− \frac{n}{2}} · ∥∇(u_{1} − u_{2})∥2^{\frac{n}{2}} )
= c · ∫_{0}^{T} ( ∥∇u_{2}∥_{2} · ∥u_{1} − u_{2}∥_{2}^{2− \frac{n}{2}} ) · ∥∇(u_{1} − u_{2})∥2^{\frac{n}{2}}
(Youngsche Ungl.) ≤ ∫_{0}^{T} ∥∇(u_{1} − u_{2})∥_{2}^{2} + c̃∫_{0}^{T} ∥∇u_{2}∥2^{\frac{4}{4−n}} · ∥u_{1} − u_{2}∥_{2}^{2}

Eingesetzt in Gleichung (60) erhält man

\frac{1}{2}∥v(T)∥_{2}^{2} + ∫_{0}^{T} ∥∇v∥_{2}^{2}

= |∫_{Ω_{t}}(v∇)u_{2} · v|
≤ ∫_{0}^{T} ∥∇v∥_{2}^{2} + c̃∫_{0}^{T} ∥∇u_{2}∥2^{\frac{4}{4−n}} · ∥u_{1} − u_{2}∥_{2}^{2}

⇒ ∥u_{1}(T) − u_{2}(T)∥_{2}^{2} ≤ c ∫_{0}^{T} ∥∇u_{2}∥2^{\frac{4}{4−n}} · ∥u_{1} − u_{2}∥_{2}^{2} =: ∫_{0}^{T} α(t)∥u_{1} − u_{2}∥_{2}^{2} (61)

Diese Ungleichung würde in allen Dimensionen gelten, falls ∥∇u_{2}∥_{∞} endlich ist, d.h. falls α(t) integrierbar ist. Im Fall n = 2 gilt für jede glatte Lösung die folgende Energie-Ungleichung

∥u(t)∥_{2}^{2} + s ∫_{0}^{1} ∥∇u(s)∥_{2}^{2} ds ≤ ∥u_{0}∥_{2}^{2} + 2 ∫_{0}^{t} ∫_{Ω} f · u dx ds ,

die an dieser Stelle die gesuchte Abschätzung leisten kann:

∫_{0}^{T} ∥∇u_{2}∥_{2}^{2} ≤ \frac{1}{2}(∥u_{2}(0)∥_{2}^{2} − ∥u_{2}(T)∥_{2}^{2}) + ∫_{Ω_{t}} f · u_{2} < ∞

Dann liefert die Gronwall-Ungleichung den letzten Schritt des Beweises:

Nach Gronwall gilt: Für integrierbare α , v mit α(T) ≥ 0 und

∥v(T)∥_{2}^{2} ≤ ∥v(0)∥_{2}^{2} + ∫_{0}^{T} α∥v∥_{2}^{2}

ist

∥v(T)∥_{2}^{2} ≤ ∥v(o)∥_{2}^{2} · e^{N(T)} ,

wobei

N(T) = ∫_{0}^{T} α(t)dt .

In unserem Fall ist v = u_{1} − u_{2} und v(0) = 0,also

⇒ ∥u_{1}(T) − u_{2}(T)∥_{2}^{2} ≤ 0 · e^{N(T)} = 0
⇒ u_{1} = u_{2} . q.e.d.

Der obige Beweis zeigt, daß die Existenz einer stetigen schwachen Lösung auf (0, T) die Existenz weiterer Lösungen auf diesem Intervall ausschließt, also zumindest lokale Eindeutigkeit.

7.2.3 Existenz und Eigenschaften starker Lösungen

Wie vorher schon erwähnt, ist es nun das Ziel, eine starke Lösung der Navier-Stokes-Gleichungen über die Integralgleichung und die Stokes-Halbgruppe zu konstruieren. Zur Vereinfachung sei im weiteren f = 0 und Ω beschränktes Gebiet im ℝ^{n} für n ≥ 3. Um die Notationen des Problems auf die allgemeine Dimension n zu erweitern, sei die Definition des Banachraumes L auf

L^{n,σ} = Vervollständigung von H in L^{n}(Ω)

erweitert.

Das folgende Theorem umfaßt alle zu zeigenden Eigenschaften einer Lösung (vergl. [39], Theorem 6.1, S. 9):

Theorem 7.1 Sei a ∈ L^{n,σ}(Ω) Anfangswert der vollen Navier-Stokes-Gleichungen. Dann gilt:

• Es existiert eine maximale Zeit T > 0 derart, daß eine Lösung u ∈ C([0, T) ; L^{n,σ}) des Anfangs-Randwert-Problems existiert, welche eindeutig und glatt ist für t > 0.

• Die Normen

sup_{t≤T_{1}} t^{\frac{1}{2}} ∥∇u∥_{n} und sup_{t≤T_{1}} t^{(1− \frac{n}{r})/2}∥u∥_{r}

sind endlich für T_{1} < T , r ≥ n.

• Ist \frac{2}{s} + \frac{n}{r} = 1 , r > n, so gilt

∫_{0}^{T_{1}} ∥u∥_{r}^{s}dt < ∞ .

• Ist die maximale Existenz-Zeit T endlich, so ist u gleichmäßig stetig auf [0, T).

• T = ∞, falls ∥a∥_{n} genügend klein ist. Die oben genannten Normen mit T_{1} = ∞ sind in diesem Falle beschränkt durch c∥a∥_{n}.

Beweis: Gemäß der bereits erwähnten Integralgleichung sei das folgende Iterationsschema definiert:

u_{j+1} =: u_{0} − F (u_{j}) ,  (62)
mit u_{0}(t) = e^{−tA}a(t) ,
F (v)(t) = ∫_{0}^{t} e^{−(t−s)A} · P (v(s)∇v(s)) ds

Wähle nun ein festes δ ∈ (0, 1] und 0 < T ≤ ∞ und definiere

K_{j} := sup_{t<T} t^{(1−δ)/2}∥u_{j}(t)∥\frac{n}{δ} ,
K_{j}^{1} := sup_{t<T} t^{1/2}∥∇u_{j}(t)∥_{n} und
R_{j} := max{K_{j} , K_{j}^{1}}

Wählt man nun p = \frac{n}{1+δ} , q = \frac{n}{δ} in den L^{p} − L^{q}-Abschätzungen (55) und (57) bezüglich der Stokes-Halbgruppe, so erhält man

K_{j+1} = sup_{t<T} t^{(1−δ)/2}∥u_{j+1}(t)∥\frac{n}{δ}
= sup_{t<T} t^{(1−δ)/2}∥u_{0} − F (u_{j})∥_{q}
= sup_{t<T} t^{(1−δ)/2}∥u_{0} − ∫_{0}^{t} e^{−(t−s)A} · P (u_{j}∇u_{j})∥_{q}
≤ sup_{t<T} t^{(1−δ)/2}∥u_{0}∥_{q} + sup_{t<T} t^{(1−δ)/2}∥ ∫_{0}^{t} e^{−(t−s)A} · P (u_{j}∇u_{j})∥_{q}
≤ K_{0} + sup_{t<T} t^{(1−δ)/2} ∫_{0}^{t} ∥e^{−(t−s)A}P (u_{j}∇u_{j})∥_{q}
(Ungl. (55)) ≤ K_{0} + sup_{t<T} t^{(1−δ)/2} ∫_{0}^{t} c_{0}(t − s)^{− \frac{n}{2} ( \frac{1+δ−δ}{n})}∥P (u_{j}∇u_{j})∥_{p}
= K_{0} + sup_{t<T} t^{(1−δ)/2}c_{0} · ∫_{0}^{t}(t − s)^{− \frac{1}{2}} ∥P (u_{j}∇u_{j})∥_{p}
(Ungl. (57)) ≤ K_{0} + sup_{t<T} t^{(1−δ)/2}c_{0} · c_{2} · ∫_{0}^{t}(t − s)^{− \frac{1}{2}} ∥u_{j}∥_{q}∥∇u_{j}∥_{n}
≤ K_{0} + c · K_{j} · K_{j}^{1}  (63)
und analog K_{j+1}^{1} ≤ K_{0}^{1} + c · K_{j} · K_{j}^{1}  (64)

Ungleichungen (64) und (65) liefern direkt

R_{j+1} = max{K_{j+1} , K_{j+1}^{1}}
≤ max{K_{0} , K_{0}^{1}} + γ · R_{j}^{2}
= R_{0} + γ · R_{j}^{2}

Setzt man nun R_{0} als genügend klein voraus, z.B. R_{0} ≤ \frac{1}{6γ}, so erhält man induktiv

R_{j} ≤ 2 R_{0} ≤ \frac{1}{3γ} .

Betrachtet man analog die Folge

w_{j}(t) := u_{j}(t) − u_{j−1}(t) , w_{0}(t) := u_{0}(t) ,

so erfüllen die analogen Werte R̃_{j} die Ungleichung

R̃_{j+1} ≤ 2γ R̃_{j} R_{j} ≤ \frac{2}{3}R̃_{j} .

Für die Folge w_{j} heißt das:

∥w_{j+1}∥_{n} = ∥u_{j+1} − u_{j}∥_{n}
= ∥ − ∫_{0}^{t} e^{−(t−s)A}P (u_{j}∇u_{j}) + ∫_{0}^{t} e^{−(t−s)A}P (u_{j−1}∇u_{j−1})∥_{n}
≤ ∫_{0}^{t} ∥e^{−(t−s)A}P (u_{j}∇u_{j}) − e^{−(t−s)A}P (u_{j−1}∇u_{j−1})∥_{n}
(Ungl. (55)) ≤ c · ∫_{0}^{t}(t − s)^{− \frac{n}{2} ( \frac{1+δ−1}{n})}∥P (u_{j}∇u_{j}) − P (u_{j−1}∇u_{j−1})∥\frac{n}{1+δ}
≤ c · ∫_{0}^{t}(t − s)^{− \frac{δ}{2}} (∥w_{j}∇u_{j}∥\frac{n}{1+δ} + ∥u_{j}∇w_{j}∥\frac{n}{1+δ})
< c · R̃_{j} · R_{j}
≤ c · (2/3)^{j}

Man sieht, daß die Folge u_{k} = ∑_{j=0}^{k} w_{j} eine Cauchyfolge in C([0, T], L^{n}(Ω)), falls die u_{j} trotz der Iterationsschritte innerhalb dieses Banachraumes verbleiben. Siehe dazu [39], Beweis zu Theorem 6.1, S. 11. In diesem Fall konvergiert die Iterationsfolge gegen eine Lösung der Integralgleichung

u(t) = e^{−tA}a − ∫_{0}^{t} e^{−(t−s)A}P (u∇u) ds .

Außerdem liefert die Ungleichung R_{j} ≤ 2R_{0} für δ = \frac{n}{r} die gesuchten gewichteten Abschätzungen:

sup_{t<T} t^{\frac{1}{2}} ∥∇u∥_{n} ≤ 2R_{0} < ∞  (65)
sup_{t<T} t^{(1− \frac{n}{r})/2}∥u∥_{r} ≤ 2R_{0} < ∞  (66)

Der Beweis, daß ∫ ∥u(t)∥_{r}^{s}dt < ∞ für \frac{2}{s} + \frac{n}{r} = 1 , n < r < ∞, benötigt einen kleinen Umweg: Man definiert eine geignete Abbildung U und einen passenden Index s, so daß das Interpolationstheorem von Marcinkiewicz (siehe z.B. [34], S. 108), das im folgenden noch näher erläutert werden soll, eine geeignete Abschätzung bezüglich der L^{s}-Norm liefern kann. Dazu seien zunächst zwei Eigenschaften einer Abbildung

U : L^{p}(ℝ^{n}) ⟶ L^{q}(ℝ)

beschrieben:

1. U heißt stark stetig genau dann, wenn es ein A > 0 gibt, so daß für alle

f ∈ L^{p}(ℝ^{n}):

∥U(f)∥_{q} ≤ A · ∥f∥_{p} .

2. U ist vom schwachen (p,q)-Typ genau dann, wenn es ein A > 0 gibt, so daß für alle f ∈ L^{p}(ℝ^{n}) und alle α > 0 gilt:

meas\{t ∈ ℝ| U(f)(t) > α\} ≤ A^{q} · (\frac{∥f∥_{p}}{α} )^{q} .

Offensichtlich folgt aus 1. direkt 2., da gilt:

(α^{q} · meas\{t ∈ ℝ| |U(f)(t)| > α\})^{\frac{1}{q}}
≤ (∫ |U(f)(t)|^{q} dt)^{\frac{1}{q}}
≤ A · ∥f∥_{p}

Diese Eigenschaften sollen nun auf die Abbildung

U(f)(t) := \{ \begin{matrix}∥e^{−tA}P f∥_{r} , & 0 ≤ t ≤ T \\ 0 , & sonst \end{matrix} \} , r > n fest

angewendet werden. Diese ist für alle f ∈ L_{p} mit p ≤ r definiert. Setzt man nun die benötigten Indizes als

p_{1} := \frac{n}{2} < p_{2} := r und
s_{1} := \frac{1}{1 − \frac{n}{2r}} < s_{2} := ∞ ,

so gilt:

U : L_{p_{2}}(ℝ^{n}) = L_{r}(ℝ^{n}) ⟶ L_{∞}(ℝ) = L_{s_{2}}(ℝ)

ist stark stetig, da |U(f)(t)| ≤ c_{2}∥f∥_{r}, denn der Operator

e^{−tA}P : L_{r} ⟶ L_{r}

ist, wie schon bekannt, beschränkt.

2.

U : L_{p_{1}}(ℝ^{n}) ⟶ L_{s_{1}}(ℝ)

ist vom schwachen (p_{1}, s_{1})-Typ, denn mittels Ungleichung (55) mit p = \frac{n}{2} = p_{1} , q = r = p_{2} folgt:

∥e^{−tA}P f∥_{r} ≤ c · t^{− \frac{n}{2} ( \frac{2}{n} − \frac{1}{r})}∥f∥\frac{n}{2}
= c · t^{−1+ \frac{n}{2r}} ∥f∥\frac{n}{2}
= c · t^{− \frac{1}{s_{1}}} ∥f∥\frac{n}{2}
⇒ U(f)(t) > α ⇔ 0 ≤ t^{1− \frac{n}{2r}} ≤ c · \frac{∥f∥\frac{n}{2}}{α}
⇔ 0 ≤ t ≤ c̃(\frac{∥f∥\frac{n}{2}}{α} )^{s_{1}} ,

so daß meas\{t | U(f)(t) > α\} ≤ c̃ · (\frac{∥f∥_{p_{1}}}{α} )^{s_{1}}.

Das Interpolationstheorem von Macinkiewicz sagt nun allgemein, daß

∥U(f)∥_{L_{s}(ℝ)} ≤ c · ∥f∥_{L_{p}(ℝ^{n})} ,

das heißt

(∫_{0}^{T} ∥e^{−tA}P f∥_{r}^{s} dt)^{\frac{1}{s}} ≤ c · ∥f∥_{L_{p}(ℝ^{n})}

mit von T unabhängiger Konstante c, falls für die Indizes p,s und den Interpolationsparameter θ gilt:

p ∈ ( p_{1} , p_{2} )  mit

\frac{1 − θ}{p_{1}} + \frac{θ}{p_{2}} = \frac{1}{p} und
\frac{1 − θ}{s_{1}} + \frac{θ}{s_{2}} = \frac{1}{s} .

Wählt man an dieser Stelle nun p = n mit n ∈ (\frac{n}{2}, r), so ergibt sich daraus folgendes θ:

\frac{1 − θ}{\frac{n}{2}} + \frac{θ}{r} = \frac{1}{n}

⇒ 2 − 2θ + \frac{n}{r}θ = 1
⇒ θ = \frac{1}{2 − \frac{n}{r}} ( < 1 ) .

Nimmt man wie gehabt s_{1} = ∞ und s_{2} = (1− \frac{n}{2r})^{−1}, so berechnet sich das zugehörige s mit

s = \frac{s_{1}}{1 − θ} = \frac{2r}{r − n}

⇒ \frac{2}{s} = \frac{r − n}{n} = 1 − \frac{n}{r}
⇒ \frac{2}{s} + \frac{n}{r} = 1 .

Wendet man nun letztendlich U speziell auf a an, so erhält man mit dem Theorem von Marcinkiewicz direkt, daß

(∫_{0}^{T} ∥ \underbrace{e^{−tA}P a}_{=u_{0}} ∥_{r}^{s} dt)^{\frac{1}{s}} ≤ c · ∥a∥_{n}

⇒ ∫_{0}^{T} ∥u_{0}(t)∥_{r}^{s}dt ≤ c̃ · ∥a∥_{n}^{s} < ∞ für \frac{2}{s} + \frac{n}{r} = 1 ,  (67)

womit der lineare Teil der gesuchten Abschätzung gezeigt ist. Den Rest liefert das Iterationsschema:

Ist r = \frac{n}{δ} , s = \frac{2}{1−δ} > 2, so gilt:

∥u_{j+1}∥\frac{n}{δ} = ∥u_{0} − F (u_{j})∥\frac{n}{δ}
= ∥u_{0} − ∫_{0}^{t} e^{−(t−s)A}P (u_{j}∇u_{j})∥\frac{n}{δ}
≤ ∥u_{0}∥\frac{n}{δ} + ∥ ∫_{0}^{t} e^{−(t−s)A}P (u_{j}∇u_{j})∥\frac{n}{δ}
(Ungl. 55) ≤ ∥u_{0}∥\frac{n}{δ} + ∫_{0}^{t} ∥e^{−(t−s)A}P (u_{j}∇u_{j})∥\frac{n}{δ}
(Ungl. 57) ≤ ∥u_{0}∥\frac{n}{δ} + c_{0} · ∫_{0}^{t}(t − s)^{− \frac{n}{2} ( \frac{1}{p} − \frac{δ}{n})}∥P (u_{j}∇u_{j})∥_{p}

(p = \frac{n}{1 + δ}) ≤ ∥u_{0}∥\frac{n}{δ} + c_{1} · ∫_{0}^{t}(t − s)^{− \frac{1}{2}} ∥u_{j}∥\frac{n}{δ} ∥∇u_{j}∥_{n}
(Ungl. 66) ≤ ∥u_{0}∥\frac{n}{δ} + c · R_{0} · ∫_{0}^{t}(t − s)^{− \frac{1}{2}} t^{− \frac{1}{2}} ∥u_{j}∥\frac{n}{δ}

Wählt man nun ein β mit 2 < β < s und wendet die Höldersche Ungleichung an, so ergibt sich daraus

∥u_{j+1}(t)∥β_{\frac{n}{δ}} ≤ c∥u_{0}(t)∥β_{\frac{n}{δ}} + cR_{0}^{β}t^{−1} ∫_{0}^{t} ∥u_{j}(τ)∥β_{\frac{n}{δ}} dτ .  (68)

Den nächsten Schritt liefert Hardy’s Ungleichung, die folgendes besagt:

Für F = ∫_{0}^{t} f dτ und p > 1 gilt

∫_{0}^{∞}(F · t^{−1})^{p}dt < (\frac{p}{p − 1})^{p} ∫_{0}^{∞} f^{p} dτ ,

so daß man aus \frac{s}{β} > 1 erhält:

(∫_{0}^{T} (\underbrace{cR_{0}^{β} ∫_{0}^{t} ∥u_{j}(τ)∥β_{\frac{n}{δ}} dτ}_{=F} ·t^{−1})^{\frac{s}{β}} )^{\frac{1}{s}}

≤ (( \frac{\frac{s}{β}}{\frac{s}{β} − 1})^{\frac{s}{β}} · ∫_{0}^{T} c̃R_{0}^{s}∥u_{j}(τ)∥s_{\frac{n}{δ}} )^{\frac{1}{s}}

= (\frac{s}{s − β} )^{β} · (c̃R_{0}^{s} ∫_{0}^{T} ∥u_{j}(τ)∥s_{\frac{n}{δ}} )^{\frac{1}{s}}

= c · R_{0} · \underbrace{(∫_{0}^{T} ∥u_{j}(τ)∥_{r}^{s}dτ)^{\frac{1}{s}}}_{=:I_{j}}

= c · R_{0} · I_{j}  (69)

Wir haben also die Ungleichung

I_{j+1} ≤ cI_{0} + cR_{0}I_{j}

für I_{j} = (∫_{0}^{T} ∥u_{j}(t)∥_{r}^{s}dt)^{\frac{1}{s}} erhalten, die uns nun durch Grenzübergang die gesuchte Abschätzung liefert.

∫_{0}^{T} ∥u(t)∥_{r}^{s}dt ≤ c ∫_{0}^{T} ∥u_{0}(t)∥_{r}^{s}dt < ∞  (70)

Also ist eine Lösung der Navier-Stokes-Gleichungen als Lösung der Integralgleichung auf dem Intervall [0,T] konstruiert worden, unter der Voraussetzung, daß R_{0} hinreichend klein ist. Hierbei war R_{0} folgendermaßen gewählt:

u_{0}(t) = e^{−tA}a(t) ,
K_{0} = sup_{t≤T} t^{(1−δ)/2}∥e^{−tA}a∥\frac{n}{δ}
(Ungl. 55) ≤ sup_{t≤T} t^{(1−δ)/2} · c_{0} · t^{− \frac{n}{2} ( \frac{1}{n} − \frac{δ}{n})}∥a∥_{n}
= c_{0} · ∥a∥_{n}
K_{0}^{1} = sup_{t≤T} t^{1/2}∥∇e^{−tA}a∥_{n}
(Ungl. 56) ≤ sup_{t≤T} t^{1/2} · c_{1} · t^{− \frac{n}{2n}} ∥a∥_{n}
= c_{1} · ∥a∥_{n}
⇒ R_{0} = max{K_{0} , K_{0}^{1}}
≤ c · ∥a∥_{n}

Da unabhängig von T gilt, daß R_{0} ≤ c∥a∥_{n}, implizieren also kleine Anfangswerte bereits die globale Existenz einer Lösung.

Andererseits genügt R_{0} auch für beliebiges T > 0 und μ > 0 (siehe [39], S. 12) der Ungleichung

R_{0} ≤ ∥a − e^{−μA}∥_{n} + ∥a∥_{n} (\frac{T}{T + μ})^{\frac{1}{2}} .

Wählt man hier zuerst ein geeignet kleines μ und dann ein passendes T, so liefert diese Ungleichung auch für große Anfangswerte eine lokale Beschränkung von R_{0} und damit die Existenz einer lokalen Lösung. Zum Schluß bleibt dann nur noch zu zeigen, daß die sogenannte milde Lösung der Integralgleichung bereits starke Lösung der Navier-Stokes-Gleichungen ist. H. Kozono und T. Ogawa zeigten in [24], S. 799, ähnlich wie M. Wiegner in [39], S. 12, daß für jede milde Lösung u der Navier-Stokes-Gleichungen gilt:

u(t) = e^{−tA}a − ∫_{0}^{t} e^{−(t−s)A}P (u∇u) ds ∈ D(A) ,
Au ∈ C_{0}((0, T) , L^{\frac{n}{1+δ}})

mit ∥Au∥_{n} ≤ ct^{−1}(log t)^{1−1/n} für genügend großes t
und u_{t} + Au + P (u∇u) = 0 ,

womit das betrachtete Existenz- und Eindeutigkeitstheorem vollendet wäre.

An dieser Stelle sind die wesentlichen theoretischen Grundlagen zum Verhalten numerisch gewonnener Lösungen der Navier-Stokes-Gleichungen geklärt. Kehrt man jetzt wieder zurück zum eigentlichen Ausgangsproblem, den künstlichen Herzklappen, so ist es also möglich, den Blutfluß durch eine künstliche Herzklappe mittels der Navier-Stokes-Gleichungen mit geeigneten Randbedingungen so exakt wie möglich zu beschreiben. Die Anwendung numerischer Lösungverfahren liefert anschließend eine eindeutige Darstellung von Geschwindigkeit und Druck an jedem Ort innerhalb des berechneten Gebietes und zu jedem betrachteten Zeitpunkt.

Der nächste Schritt ist jetzt die Implementierung der Gleichungen, der passenden Randbedingungen, einer geeigneten Geometrie des zu diskretisierenden Gebietes und eines Finite-Elemente-Lösers, der in der Lage ist, das nötige Gleichungssystem aufzustellen und effektiv zu lösen.

8 Implementierung unter “Diffpack 3.0”

Das folgende Kapitel soll eine Möglichkeit zur Implementierung eines Finite-Elemente-Lösers (FE-Lösers) für die Navier-Stokes Gleichungen aufzeigen. Für den Aufbau eines geeigneten FE-Lösers wurde “Diffpack 3.0” benutzt, eine flexible und vielfältige Software, die speziell zur Lösung partieller Differentialgleichungen mit der Methode der finiten Elemente konstruiert wurde.

Diffpack 3.0 ist eine Sammlung von Bibliotheken, die alle notwendigen Bausteine zur Lösung komplexer Probleme und Aufgaben enthalten: Lineare und nichtlineare Gleichungssysteme, numerische Lösungsmethoden für diese Systeme, Gittererzeuger, skalare Felder und Vektorfelder auf Gittern, verschiedene Elementtypen und letztendlich einige Tools zur Visualisierung der Ergebnisse. Diffpack 3.0 ist in C++ geschrieben, was die Vor- und Nachteile der objektorientierten Programmierung mit sich bringt. Die Details werden bei genauerer Betrachtung des Programm-Codes deutlich.

Das Ziel dieses Kapitels ist es, eine Lösung der stationären, inkompressiblen, 2-dimensionalen Navier-Stokes Gleichungen auf der Geometrie einer mechanischen Aortenklappe zu finden und zu visualisieren. Die zu lösenden Gleichungen sind also

(u · ∇) u = −grad p + \frac{1}{Re}Δu auf dem ganzen Gebiet Ω ,
div u = 0 in ganz Ω und
u = 0 auf dem Rand ∂Ω .

8.1 Geometrie

Das Gebiet an einer mechanischen, zweiflügeligen Herzklappe weißt mehrere Symmetrieebenen auf, wie z. B. aus Abbildung 12, S. 26, entnehmbar ist. Ebenso wie in der Doktorarbeit von M. J. King, die im folgenden als Leitfaden zur Modellierung dient, sollen also nur auf der oberen Hälfte der Klappe und des folgenden Aortenabschnittes Geschwindigkeit und Druck berechnet werden. Auf diesem Gebiet soll bei unterschiedlichen Klappenstellungen bzw. variierenden Öffnungswinkeln (75^{∘}, 80^{∘}, 85^{∘} und 90^{∘}) ein Finite-Elemente Gitter erzeugt werden. Die Penalty-Methode, die sowohl hier als auch in der Arbeit von King verwendet wurde und im folgenden noch näher erläutert wird, liefert beste Ergebnisse auf quadratischen oder rechteckigen Elementen und unakzeptable Werte auf Triangulierungen. Viele gängige Gitter-Erzeugungs-Tools sind damit für diese Zwecke nicht geeignet und man muß z. B.

Abbildung 15: Abstrahiertes, 2-dimensionales Gebiet zur Beschreibung einer zweiflügeligen, mechanischen Herzklappe. Die weitere Modellierung beschränkt sich auf die obere Hälfte; Größenangaben aus [23], S. 44ff.

auf die Diffpack-internen Methoden zurück greifen. Das stark abstrahierte, zur Modellierung verwendete Gebiet ist in Abbildung 15 skizziert.

Zur Erzeugung eines geeigneten Gitters wurden zwei unterschiedliche Methoden verwendet: Einerseits findet man unter den Diffpack-Beispielen eine Geometrie mit dem bezeichnenden Namen “BOX WITH BELL”, die ein 2-dimensionales Rohr mit einer Ausbuchtung in Form einer Gauss-Kurve beinhaltet. Benutzt man hier die unter Diffpack bereits implementierte Funktion “addBoIndsNodes”, zur Definition zusätzlicher Randpunkte mit den dazugehörigen Randindikatoren, so kann man das in Abbildung 15 gezeigte Gebiet recht einfach konstruieren. (siehe Anhang B.6) Der Nachteil dieser Methode besteht darin, daß die Funktion “addBoIndsNodes” als Eingabe mit Intervallen in x- und in y-Richtung arbeitet, so daß eine Schrägstellung des Klappenflügels innerhalb dieser Methode nicht möglich war.

Um komplexere Gitter zu generieren, kann man die in Diffpack bestehende Super-Elemente-Methode benutzen, die das Gebiet in geometrisch einfache Teilgebiete gliedert, diese jeweils mit einer adäquaten Zerlegung versieht und schließlich ein komplettes Gitter zusammensetzt. (siehe dazu [27], Kapitel “Super-Elements”) Diese Methode hat einige sehr empfindliche Parameter, die man besonders beachten sollte, wenn man den Übergang zwischen den einzelnen Super-Elementen gestaltet:

Übergibt man einem Super-Element die Randdaten des Nachbar-Elementes, so kann man zu viele oder auch zu wenige Angaben machen, was direkt zur Divergenz des iterativen Lösers führen kann. Leider sind die Gesetztmäßigkeiten hier nicht offensichtlich, so daß auch die in dieser Arbeit gelieferten Gitter und daraus resultierenden Werte nicht optimal sind. Aber die Methode der Super-Elemente birgt eine große Ausbaufähigkeit und Flexibilität in sich. Auf diese Weise sind z. B. auch verschiedene Öffnungswinkel der Klappenflügel realisierbar. Das Gitter mit Öffnungswinkel 75^{∘} hat folgende Gestalt:

Abbildung 16: FE-Gitter bei Öffnugswinkel 75^{∘} unter Verwendung der Super-Elemente-Methode

Zu weiteren Details zur Gittergenerierung siehe Anhang B.4.

8.2 Penalty-Methode

Ein Hauptproblem bei der Lösung der Navier-Stokes Gleichungen ist die Berechnung der Druckes p bzw. die Behandlung des Termes grad p dar. Eine zwar sehr empfindliche, aber zur Implementierung gut geeignete Methode ist die Penalty-Methode, die die Kontinuitätsgleichung approximiert durch:

div u = −λ^{−1} · p für λ → ∞ , λ ∈ ℝ
⇔ p = −λ · div u ,

so daß der Druckterm aus den zu lösenden Gleichungen eliminiert werden kann. Diese auf den ersten Blick sehr einfach erscheinende Methode erfordert einige Vorsichtsmaßnahmen: Einerseits sollte für die λ-Terme in den folgenden Galerkin-Schritten eine andere Integrationsregel verwendet werden, als für die λ-freien Terme (siehe dazu [26], S. 431, und Programmcode) und andererseits ist die Konditionszahl der Koeffizientenmatrix im zu lösenden linearen System proportional zu λ. Wählt man also ein großes λ zur guten Approximation der Gleichung div u = 0, so sinkt die Konvergenzgeschwindigkeit der iterativen Löser stark. Da der Effekt der üblichen Vorkonditionierungsmethoden auch nur verschwindend gering ist, wird das lineare Gleichungssystem in jedem nichtlinearen Newton-Raphson- oder Successive-Substitution-Schritt durch direkte Methoden wie die Gauss-Elimination gelöst. Durch diese Schwierigkeiten sind Simulations-Programme, die auf der Penalty-Methode basieren, auf Gebiete bzw. Gitter mit wenigen Tausend Knotenpunkten beschränkt.

Setzt man also den approximierten λ-Term in die Gleichungen ein, so erhält man:

(u · ∇) u = λ grad(div u) + \frac{1}{Re}Δ u und
div u = −λ^{−1} · p .

Verwendet man nun statt der Reynolds-Zahl explizit die Dichte ρ und die Viskosität μ, so erhält man für ein geeignet gewähltes λ:

ρ (u · ∇) u − μ Δ u − λ grad(div u) = 0

Die Methode der gewichteten Residuen bzw. hier die Galerkin-Methode (siehe Kapitel 6) liefert nun durch Integration über das ganze Gebiet Ω die zu implementierende dikretisierte Gleichung. Wie schon zuvor seien die Basisfunktionen mit N_{i} bezeichnet, so daß die Approximation der Geschwindigkeit u die Gestalt

ũ_{k} = ∑_{j=1}^{n} u_{j}^{k} N_{j}

erhält, wobei die (u_{j}^{k}) also die gesuchten Koeffizienten für j = 1, …, n und k = 1, 2 (bei 2 räumlichen Dimensionen) sind. Mit den gleichen N_{i} als Wichtungsfunktionen, unter Berücksichtigung der Randbedingungen und mittels partieller Integration ensteht das folgende Gleichungssystem:

(∫_{Ω} ρ · N_{i} · ũ_{k} · ∇N_{j} + ∫_{Ω} μ · ∇N_{i} · ∇N_{j} + ∫_{Ω} λ · ∇N_{i} · ∇ũ_{k}) · (u_{j}^{k}) = 0

für i, j = 1, …, n und k = 1, 2.

8.3 Lösung des nichtlinearen Systemes

Im folgenden sollen die beiden verwendeten Methoden zur Lösung des oben beschriebenen nichtlinearen Systemes genauer erläutert werden. Beide werden getrennt voneinander auf die λ-Terme (siehe im Programmcode Funktion intergrandsNonReduced) und die λ-freien Terme (siehe Funktion integrandsReduced) angewendet.

8.3.1 Successive Substitutions

Eine Methode zur Lösung nichtlinearer Gleichungssysteme wie das oben beschriebene ist die einfach Iterationsmethode der “Successive Substitutions”, auch bekannt unter dem Namen “Picard Iterationen”. Anstelle des nichtlinearen Systemes der Form

A(v) · v = b

wird in jedem Iterationsschritt das System

A(v^{k}) · v^{k+1} = b

für k = 0, 1, 2, … gelöst, bis der Fehler ∥u^{k+1}−u^{k}∥ genügend klein ist. Ein geeigneter Anfangsvektor ist meist v^{0} = 0. Falls dieses Methode nicht konvergiert, so hilft in den meisten Fällen ein Ralaxationsparameter ω ∈ (0, 1):

Wie schon bei anderen iterativen Verfahren beschrieben, wird zuerst eine neue Approximation v^{∗} aus A(v^{k}) · v^{∗} = b berechnet. Den neuen Wert v^{k+1} gewinnt man dann aus dem gewichteten Mittel von v^{k} und v^{∗}:

v^{k+1} = ω v^{∗} + (1 − ω) v^{k} .

8.3.2 Newton-Raphson-Methode

Zur Erläuterung der “Newton-Raphson-Methode” betrachten wir zunächst nur eine nichtlineare Gleichung F (v) = 0 mit nur einer Unbekannten v. Die zugrunde liegende Idee ist, das Gleichungssystem F (v) in Abhängigkeit von einer leichter zu berechnenden approximativen Lösung v^{k} anzunähern. Man löst also ein lieneares Gleichungssystem M(v; v^{k}) ≈ F (v) mit der Eigenschaft, daß M(v; v^{k}) = 0 einfach zu lösen ist. Die Form dieses linearen Gleichungssytemes erinnert an den linearen Teil der Taylor-Entwicklung von F im Punkt v = v^{k}:

M(v; v^{k}) = F (v^{k}) + \frac{dF}{dv} (v^{k})(v − v^{k}) .

Ist die nächste iterativ gewonnene Lösung v^{k+1} nun die Lösung von M(v; v^{k}) = 0, so gilt also

v^{k+1} = v^{k} − \frac{F (v^{k})}{\frac{dF}{dv} (v^{k})} .

Dieses Schema kann nun recht einfach auf mehrdimensionale Probleme übertragen werden. Wieder approximiert man F (v) durch eine lineare Funktion M(v; v^{k}) in Abhängigkeit von einer bestehenden Approximation v^{k} zu v:

M(v; v^{k}) = F (v^{k}) + J(v^{k}) · (v − v^{k}) ,

wobei J = ∇F die Jacobi-Matrix von F bezeichnet. Zur Verbesserung der Konvergenz kann wie bei den “Successive Substitutions” ein Relaxationsparameter ω ∈ (0, 1) eingefügt werden.

8.4 Programmaufbau

Zur Lösung des betrachteten Problemes wurde ein unter Diffpack bereits vorhandener FE-Löser der Navier-Stokes Gleichungen auf das betrachtete Problem angepaßt. Da nicht alle darin verwendeten Klassen für den Benutzer offensichtlich und frei zugänglich sind, war es in der vorhandenen Zeit nicht möglich, die in diesem Fall überflüssige Zeitschleife zu eliminieren.

Das Simulationsprogramm NsPenalty.cpp bzw. der Programmcode hat folgenden Aufbau:

  → define

Definition der Items des GUI-Menüs

  → scan

Einlesung der Parameter aus dem GUI-Menü: Gitter, Randindikatoren, Viskosität, Dichte, Penalty-Parameter usw. ⇒ Initialisierung der Variablen

  → fillEssBC

Randindikatoren:

bo-ind 1: inlet-Geschwindigkeitsfeld

bo-ind 2: u=0

bo-ind 3: v=0

bo-ind 4: w=0

⇒ Zuordnung der Randwerte zu den entsprechenden Randpunkten

  → solveProblem

Aufruf des eigentlichen Lösers mit innerer Zeitschleife

      → setIC

Anfangsbedingungen werden auf 0 gesetzt.

      → timeLoop

Zeitschleife: In jedem nächsten Zeitschritt wird Funktion solveAtThisTimeStep aufgerufen!

→ inletVelocity

Randbedingung bei x=0: parabolisches Inlet-Geschwindigkeits-Profil

→ integrands

Berechnung des Druckes p mittels der Penalty-Methode

→ solveAtThisTimeStep

Aufruf der Funktion makeAndSolveLinearSystem mit passenden Randbedingungen zum gegenwärtigen Zeitpunkt

→ makeAndSolveLinearSystem

Initialisierung und Lösung des eigentlichen globalen Gleichungssystemes

→ calcElmMatVec

Berechnung und Zusammensetztung der Element-Matrizen und -Vektoren

→ integrands

Berechnung der einzelnen Einträge für alle Elemente mit jeweiliger Integrationsregel

→ integrandsNonReduced

Berechnung der λ-freien Terme, Implementierung von Newton Raphson und Successive Substitution

→ intergrandsReduced

Berechnung der λ-Terme mit Newton Raphson oder Successive Substitution

Details können dem Programmcode dem Anhang B entnommen werden. Zu weiteren Informationen zu den unter Diffpack benutzten C++-Klassen siehe [26], S. 431ff.

8.5 Ergebnisse

Die Lösungen für Druck p bzw. Geschwindigkeit v bei unterschiedlichen Öffnungswinkeln der Klappenflügel werden als Farb-Plot bzw. als Vektor-Plot ausgegeben. Die folgenden Abbildungen sind mit dem in Diffpack zu findenden Visualisierungs-Too “Plotmtv” entstanden. Siehe dazu [26], S. 232ff. Es wird stets Blutfluß von links nach rechts betrachtet.

Abbildung 17: Druck p auf BOX WTH BELL-Geometrie bei Öffnungswinkel 90^{∘}

Betrachtet man die visualisierten Lösungen des FE-Lösers für p und v, so kann man deutlich die Wirbelbildungen in der Sinus-Ausbuchtung und in der Nähe der Klappenflügel erkennen. Verfolgt man nun die Entwicklung der Druck- und Geschwindigkeits-Ergebnisse bei abnehmendem Öffnungswinkel von 90^{∘} bis 75^{∘}, dann wird deutlich, daß sich außer im Sinus noch Wirbelbildungen am unteren, linken Rand des Klappenflügels finden lassen, durch die eine Region mit leichtem Unterdruck entsteht. Dieser trägt im Realfall maßgebend zur Bewegung des Klappenflügels bei, indem der “Sog” am linken Ende und der Druck durch den Sinus-Wirbel am rechten Ende des Klappenflügels diesen zurück schwingen lassen. Bei einem Öffnungswinkel von 90^{∘} ist diese Wirbelbildung direkt an dem Klappenflügel nicht zu sehen, während sie bei 75^{∘} stark ausgeprägt ist.

Im Vergleich mit den Ergebnissen von Mary J. King stellt man fest, daß auch dieses stark vereinfachte Modell schon erahnen läßt, was King mit einem 3-dimensionalen, zeitabhängigen Modell gezeigt hat: Das ideale Verhältnis zwischen den unterschiedlichen Kräften, die auf die Klappenflügel einwirken, ist bei einem Öffnungswinkel Abbildung 18: Geschwindigkeit v auf BOX WTH BELL-Geometrie bei Öffnungswinkel 90^{∘}

Abbildung 20: Geschwindigkeit v auf Superelemente-Geometrie bei Öffnungswinkel 90^{∘}

Abbildung 19: Druck p auf Superelemente-Geometrie bei Öffnungswinkel 90^{∘} Abbildung 21: Druck p auf Superelemente-Geometrie bei Öffnungswinkel 85^{∘} Abbildung 23: Druck p auf Superelemente-Geometrie bei Öffnungswinkel 80^{∘} Abbildung 25: Druck p auf Superelemente-Geometrie bei Öffnungswinkel 75^{∘} von ca. 85^{∘} erreicht, was der Konstruktion einer St. Jude-Herzklappenprothese entspricht. Ein maximaler Öffnungswinkel von mehr als 85^{∘} kann das reibungslose Schließen der Klappe gefährden, während ein maximaler Öffnungswinkel von weniger als 85^{∘} zu starke Turbulenzen verursacht, die sowohl dem Klappenmaterial als auch dem umliegenden Gewebe schaden, wie schon in Kapitel 2.2 beschrieben.

Zusammenfassung

Der Überblick über die biologischen, physikalischen und mathematischen Gegebenheiten zeigt, wie viele Möglichkeiten und ungeklärte Fragen man auf dem Gebiet der numerischen Simulationen noch antreffen kann. Modelliert man den Fluß einer viskosen Flüssigkeit wie den des menschlichen Blutes, so stößt man unweigerlich auf die Navier-Stokes Gleichungen, die mathematisch einige Schwierigkeiten beinhalten. Da es nicht möglich ist, eine analytische Lösung zu berechnen, ist man auf numerische Verfahren wie die Finite-Elemente-Methode zusammen mit dem Galerkin-Verfahren angewiesen. Das Hauptproblem der Navier-Stokes Gleichungen ist also ein Beweis der Existenz und Eindeutigkeit einer Lösung in mindestens drei Dimensionen, um den reell betrachteten Problemstellungen gerecht zu werden.

An die Klärung der mathematischen Hintergründe kann sich dann die Implementierung eines geeigneten Finite-Elemente-Lösers anschließen, wie es zum Beispiel mittels der C++-Bibliothek “Diffpack 3.0” möglich ist. In dieser Arbeit ist es aus Zeitgründen nicht vollständig gelungen, den vorhandenen Löser zu optimieren, d. h. es treten je nach Wahl der möglichen Parameter Konvergenzprobleme bei den iterativen Lösern des nichtlinearen globalen Gleichungssystemes, welches aus der Methode der Finite-Elemente resultiert, auf. Ein entscheidender Parameter hierbei ist die Reynoldszahl, die den Effekt von Viskosität auf das Flußverhalten beschreibt. Zur Behebung dieses Problemes könnte man für die Reynoldszahl einen “Continuation”-Parameter (siehe [26], S. 351) in die Gleichungen einfügen. Das bedeutet, man berechnet zunächst eine Lösung bezüglich einer sehr kleinen, physiologisch nicht angemessenen Reynoldszahl, um dieses Lösung als Startvektor für den nichtlinearen Löser in einer weiteren iterativen Schleife zu benutzen. Im Laufe der Iterationsschritte wird die Reynoldszahl dann durch den “Continuation”-Parameter auf die gewünschte Größe angehoben.

Trotz der beschriebenen Schwierigkeiten und starken Vereinfachungen, läßt sich der Blutfluß durch eine zweiflügelige, mechanische Herzklappe mit dem angegebenen Simulationsprogramm berechnen und skizzieren. Die Simulation des Flußverhaltens bzw. die Betrachtung von Geschwindigkeit und Druck ermöglicht einen detaillierten Einblick in die Strömungsverhältnisse, der Aufschluß über Schwachstellen in der Konstruktion einer solchen Klappenprothese geben kann. Schon kleine Änderungen der Maße von Klappenring, Klappenflügel und besonders des maximalen Öffnungswinkels führen zu einem anderen Flußbild, welches andere Auswirkungen haben kann. Dazu zählen Nebeneffekte wie z. B. Wucherungen der Aortenwand durch zu starken Aufprall des Blutstromes oder Erosionerscheinungen an den Klappenflügeln durch zu hohe Druckunterschiede. Mit der Hilfe eines Simulationsprogrammes können einzelne Parameter verändert und die Folgen direkt beobachtet werden.

Ein maximaler Öffnungswinkel von ca. 85^{∘} scheint im betrachteten Fall optimal zu sein.

Zur genauen Lokalisierung von Wirbeln und Turbulenzen ist allerdings ein 3-dimen-sionales, zeitabhängiges Modell nötig, das besser auf die sich ändernden Randbedingungen des Herzzyklus eingehen kann. Ergebnisse eines solchen Simulations-Programmes sind z. B. in der Arbeit von Mary J. King (s. [23]) zu finden. Der nächste Schritt ist nun ein Modell, welches den Blutfluß und die Bewegung der Klappenflügel als gekoppelte Probleme betrachtet, so daß man von tatsächlicher Interaktion zwischen Fluid, starrem Prothesen-Material und elastischer Außenwand sprechen kann. Doch selbst dann ist dieses Thema in naher Zukunft nicht ausgeschöpft, da die exakte Modellierung des Herzzyklus und des umliegenden Gewebes noch in weiter Ferne liegen.

A Medizinische Fachbegriffe

Anaemie sog. Blutarmut; Verminderung des Anteils der roten Blutkörperchen im Vollblut unter den Normalwert

Aorta die von der linken Herzhälfte abgehende Hauptschlagader als Stammgefäß des Körperkreislaufes, des sogenannten ”großen“ Blutkreislaufes; mit elastisch (durch entsprechenden Wandaufbau) bedingter Windkesselfunktion

Atrium Herzvorhof

Bikuspidalität die - seltene - Fehlbildung der Aorten- oder Pulmonalklappe als nur aus 2 Klappenflügeln bestehende Gebilde

distal anatomische Richtungsweisung; weiter entfernt von der Körpermitte, Herzmitte o.ä.; Gegensatz: proximal

Embolie plötzlicher Verschluß eines Blutgefäßes (meist Arterie) durch von der Blutbahn verschlepptes Gewebe, Fremdkörper oder Luft

Endokard die alle Hohlräume (einschließlich der Herzklappen) auskleidende glatte Innenhaut des Herzens, schichtweise aufgebaut - von innen zum Herzmuskel hin - aus Endothel, feinfaserigem kollagenem Bindegewebe, elastischen Fasern und einzelnen glatten Muskelzellen

Endokarditis Entzündung der Herzinnenhaut; meist als Entzündung der Herzklappen, genauer am Schließungsrand einer Klappe, aber auch im Bereich der Vorhof- und Kammerwände, Sehnenfäden und Papillarmuskeln

Epikard das mit der äußeren Oberfläche des Herzmuskels verwachsene ”innere Blatt“ des Herzbeutels als äußerste Schicht der Herzwand; erstreckt sich auch über die herznahen Teile der großen Gefäße

Haemolyse die Auflösung (Zerstörung) der roten Blutkörperchen, z.B. in Folge einer zu hohen mechanischen Beanspruchung an Herzklappenprothesen

Herzinsuffiziens akutes oder chronisches Unvermögen des Herzens, bei Belastung (=Belastungsinsuffiziens) oder schon in Ruhe (=Ruheinsuffiziens) den für den Stoffwechsel erforderlichen Blutauswurf aufzubringen

Herzminutenvolumen HMV, das vom Herzen pro Minute ausgeworfene Blutvolumen

Hypertrophie Größenzunahme eines Gewebes oder Organs nur durch Zellvergrößerung bzw. -wucherung bei normal bleibender Zellzahl intrakardial in einer oder in eine Herzhöhle Myocard der Herzmuskel; die mittlere, zwischen Endo- und Epikard gelegene Her-

zwandschicht Perikard Herzbeutel: Schutz- und Gleithülle des Herzens Rheumatisches Fieber spezifische Entzündungsreaktion auf Toxine der Strepto-

kokken serös das Blutserum betreffend; hier: serumhaltig, z.B. seröse Höhlen, wie die Bauchhöhle Sinus latein.: Bucht, Tasche; taschenartige Körperhöhle oder Organausbuchtung;

z.B. Sinus aortae die nahezu intrakardiale Ausweitung zwischen jedem der drei Aortenklappenflügeln und der Aortenwand

Stenose angeborene oder erworbene dauerhafte Einengung eines Kanals, z.B. Herz-

klappenstenose Stenosegeräusch Herzgeräusche (evt. mit tastbarem Schwirren) bei Herzklappen-

stenosen Thorax Brustkorb Thromboembolie Embolie durch Verschleppung von Gerinselteilen, z.B. bei Her-

klappenerkrankungen Thrombose die intravitale ”Blutpfropfbildung“ im Kreilaufsystem durch Aggrega-

tion von Thrombozyten Thrombozyten Blutplättchen; vom Knochenmark abstammende, kernlose Blut-

elemente Truncus pulmonalis von der rechten Herzhälfte abgehendes Stammgefäß des Lun-

genkreislaufes, des sogenannten ”kleinen“ Blutkreislaufes Ventrikel Herzkammer

B Programm-Code und Input-Files

B.1 NsPenalty1.h
#ifndef NsPenalty1_h_IS_INCLUDED
#define NsPenalty1_h_IS_INCLUDED
#include <FEM.h>                // finite element toolbox
#include <DegFreeFE.h>          // field dof <-> linear system dof
#include <LinEqAdmFE.h>         // linear solver toolbox
#include <TimePrm.h>            // time discretization parameters
#include <NonLinEqSolver.h>     // nonlinear solver interface
#include <NonLinEqSolver_prm.h> // parameters for nonlinear solvers
#include <NonLinEqSolverUDC.h>  // user's definition of nonlinear PDEs
#include <SimCase.h>            // user's definition of menu items
#include <SaveSimRes.h>         // used in storing and plotting of
                               // results

class NsPenalty1: public FEM, public NonLinEqSolverUDC
{
public:
  Handle(GridFE)         grid;    // Finite-Elemente-Gitter
  Handle(DegFreeFE)      dof;     // Zeiger: matrix dof <-> u dof
  Handle(FieldsFE)       u;       // Geschwindigkeitsfeld
  Handle(FieldsFE)       u_prev;  // u auf vorherigem Zeitlevel
  Handle(FieldFE)        p;       // Druckfeld
  Handle(TimePrm)        tip;     // Parameter zur zeitlichen
                                 // Diskretisierung
  int  inlet_profile;    // Flußprofil bei x=0
  real theta;            // "theta"-Regel zur Zeit-Diskretisierung
  real mu;               // Viskosität
  real density;          // Dichte
  real inlet_velocity;   // Flußgeschwindigkeit bei x=0
  real lambda;           // Parameter für Penalty-Methode
  Vec(real)                  nonlin_solution; // Lösung des nichtlinearen
                                             // Systemes
  Vec(real)                  linear_solution; // Lösung des linearen

                                           // Unter-Systemes
Handle(NonLinEqSolver_prm) nlsolver_prm;    // Parameter für nlsolver
Handle(NonLinEqSolver)     nlsolver;        // nichtlineare Löser
Handle(LinEqAdmFE)         lineq;           // lineare Löser
Handle(SaveSimRes)         database;        // berechnete Daten
// Zur Vermeidung von wiederholter Einlesung in die
// Funktion integrands:
Ptv(real) u_pt;                  // u am aktuellen Integrationspunkt
Ptv(real) up_pt;                 // u_prev am aktuellen
                                // Integrationspunkt
VecSimple(Ptv(real)) gradu_pt;   // grad u am aktuellen
                                // Integrationspkt.
VecSimple(Ptv(real)) gradup_pt;  // grad u_prev am aktuellen
                                // Integrationspkt.
NsPenalty1 () {};
~NsPenalty1 () {};
virtual void define (MenuSystem& menu, int level = MAIN);
virtual void scan   ();
virtual void adm    (MenuSystem& menu);
virtual void solveProblem ();
virtual void resultReport ();
virtual void saveResults ();
virtual void calcDerivedQuantities ();
virtual void fillEssBC ();
protected:
virtual void inletVelocity (Ptv(real)& v, const Ptv(real)& x);
virtual void setIC ();
virtual void timeLoop ();
virtual void solveAtThisTimeStep ();
virtual void makeAndSolveLinearSystem ();
enum Integrand_type { LAMBDA_TERMS, ORDINARY_TERMS };
Integrand_type integrands_tp;
virtual void calcElmMatVec
  (int elm_no, ElmMatVec& elmat, FiniteElement& fe);

  virtual void integrands
          (ElmMatVec& elmat, const FiniteElement& fe);
  virtual void integrandsReduced
          (ElmMatVec& elmat, const FiniteElement& fe);
  virtual void integrandsNonReduced
          (ElmMatVec& elmat, const FiniteElement& fe);
  friend class PressureIntg;
};

class PressureIntg: public IntegrandCalc
{
  NsPenalty1* data; // Zugang zu Input-Daten und Ergebnissen
public:
  PressureIntg (NsPenalty1* data_) { data = data_; }
  virtual void integrands(ElmMatVec& elmat,const FiniteElement& fe);
};
#endif

B.2 NsPenalty.cpp
#include <NsPenalty1.h>
#include <ElmMatVec.h>
#include <FiniteElement.h>
#include <readOrMakeGrid.h>
#include <DegFreeFE.h>
#include <VecSimple_Ptv_real.h>
#include <MatDiag_real.h>
#include <Puttonen.h>
void NsPenalty1:: define (MenuSystem& menu, int level)
{
  menu.addItem (level, "casename", " ", "test1");
  menu.addItem (level, "gridfile",
                "filename", "gitter/test1.grid");
  menu.addItem (level, "redefine boundary indicators",
                "GridFE::redefineBoInds syntax",
                "nb=2 names D1 d2 1=(1 3) 2=(2 4)");
  menu.addItem (level, "add boundary nodes",
                "GridFE::addBoIndNodes syntax",

              "n=1 b1=[0,0]x[0,0]");
menu.addItem (level, "time integration parameters",
              "TimePrm::scan syntax", "dt=1");
menu.addItem (level, "inlet profile", " ", "1");
menu.addItem (level, "characteristic inlet velocity",
              "used in inlet profile expressions", "25.0e-2");
menu.addItem (level, "viscosity", "mu", "4.0e-3");
menu.addItem (level, "density", "rho", "1000");
menu.addItem (level, "penalty parameter", "lambda", "10000");
menu.addItem (level, "bandwidth reduction", "for GaussElim", "ON");
LinEqAdmFE        ::defineStatic (menu, level+1);
NonLinEqSolver_prm::defineStatic (menu, level+1);
SaveSimRes        ::defineStatic (menu, level+1);
FEM               ::defineStatic (menu, level+1);
}
void NsPenalty1:: scan ()
{
MenuSystem& menu = *SimCase::menu_system;
String gridfile = menu.get ("gridfile");
grid.rebind (new GridFE ());
readOrMakeGrid (*grid, gridfile);
                  // Gittergenerierung
String redef = menu.get ("redefine boundary indicators");
 if (!redef.contains ("NONE"))
  grid->redefineBoInds (redef);
                 // neue Randwertdefinierung
String add = menu.get ("add boundary nodes");
 if (!add.contains ("NONE"))
  grid->addBoIndNodes (add);
                // zusätzliche Randpunkte
tip.rebind (new TimePrm());
tip->scan (menu.get ("time integration parameters"));
mu = menu.get ("viscosity").getReal();
lambda = menu.get ("penalty parameter").getReal();
density = menu.get ("density").getReal();

const nsd = grid->getNoSpaceDim();
database.rebind(new SaveSimRes);
database->scan (menu, nsd);
inlet_profile  = menu.get ("inlet profile").getInt();
inlet_velocity = menu.get ("characteristic inlet
                           velocity").getReal();
bool bandwidth_reduction = menu.get ("bandwidth
                                   reduction").getBool();
if (bandwidth_reduction) {
  s_o << "\nReducing the bandwidth by
           using class Puttonen.\n";
  Puttonen bwr;  bwr.renumberNodes (*grid);
 }
lineq.rebind(new LinEqAdmFE);
lineq->scan (menu);

  // Initialisierung der Variablen
p.rebind (new FieldFE (*grid,"p"));
u.rebind (new FieldsFE (*grid,"v"));
u_prev.rebind (new FieldsFE (*grid, "v_prev"));
dof.rebind (new DegFreeFE (*grid, nsd));
const int lineq_vec_length = u->getNoValues();
nonlin_solution.redim (lineq_vec_length);
linear_solution.redim (lineq_vec_length);
lineq->attach (linear_solution);
nlsolver_prm.rebind (NonLinEqSolver_prm::construct());
nlsolver_prm->scan(menu);
nlsolver.rebind (nlsolver_prm->create());
nlsolver->attachUserCode (*this);
nlsolver->attachNonLinSol (nonlin_solution);
nlsolver->attachLinSol (linear_solution);
u_pt.    redim (nsd);  up_pt.    redim (nsd);
gradu_pt.redim (nsd);  gradup_pt.redim (nsd);
}

void NsPenalty1:: adm (MenuSystem& menu)
{

SimCase::attach(menu);
define (menu); menu.prompt (); scan ();
}

void NsPenalty1:: fillEssBC ()
{
/*   Randindikatoren:
 bo-ind 1: inlet-Geschwindigkeitsfeld
 bo-ind 2: outlet-Geschwindigkeitsfeld mit
           Normalenableitungen = 0
 bo-ind 3: u=0
 bo-ind 4: v=0
 bo-ind 5: w=0
*/
dof->initEssBC ();
const int nno = grid->getNoNodes();
const int nsd = grid->getNoSpaceDim();
Ptv(real) v (nsd);
Ptv(real) x (nsd);
int i, k;
for (i = 1; i <= nno; i++) {
  if (grid->boNode (i,1)) {
    // Aufruf des inlet-Profiles
    grid->getCoor (x,i);
    inletVelocity (v, x);
    for (k=1; k <= nsd; k++)
      // k-te Komponente des i-ten
      // Knotenpunktes bekommt Wert v(k):
      dof->fillEssBC (i,k, v(k));
  }
  for (k=1; k <= nsd; k++) {
    if (grid->boNode (i, 2+k))
      dof->fillEssBC (i,k, 0.0);
  }
}
}
void NsPenalty1:: solveProblem ()

{
s_o->setRealFormat ("%13.6e");
timeLoop ();
}
void NsPenalty1:: setIC ()
{
const int nno = grid->getNoNodes();
const int nsd = grid->getNoSpaceDim();
int i,k;
for (i=1; i <= nno; i++){
  for (k=1; k <= nsd; k++) {
     u_prev()(k).values()(i) = 0.0;
     u()(k).values()(i) = 0.0;
   }
}
}
void NsPenalty1:: timeLoop ()
{
tip->initTimeLoop();
setIC ();
while (!tip->finished())
  {
    tip->increaseTime();
    solveAtThisTimeStep ();
    calcDerivedQuantities ();
    saveResults ();
    *u_prev = *u;
  }
}

void NsPenalty1:: inletVelocity (Ptv(real)& v, const Ptv(real)& x)
{
if (inlet_profile == 1)
  {
   // parabolisches Inlet-Profil in

     // x=0 mit Max. bei y=0 und Min.  bei y=0.02
      v = 0.0;
      v(1) = inlet_velocity * (1 - sqr(x(2)/0.02));
    }
 else
     errorFP("NsPenalty1::inletVelocity",
             "inlet_velocity=%d not impl.",
             inlet_velocity);
}
void NsPenalty1:: saveResults ()
{
  database->dump (*p, tip.getPtr());
  database->dump (*u, tip.getPtr());
}
void NsPenalty1:: resultReport ()
{
  // leere Funktion
}
void PressureIntg:: integrands
  (ElmMatVec& elmat, const FiniteElement& fe)
{
  const int nsd = fe.getNoSpaceDim();
  real div_v=0;
  for (int k = 1; k <= nsd; k++) {
    data->u()(k).derivativeFEM (data->gradu_pt(k), fe);
    div_v += data->gradu_pt(k)(k);
   }
  const real pressure = - data->lambda*div_v;
  const int nbf = fe.getNoBasisFunc();
  const real detJxW = fe.detJxW();
  for (int i = 1; i <= nbf; i++)
    elmat.b(i) += pressure*fe.N(i)*detJxW;
}
void NsPenalty1:: calcDerivedQuantities ()
{
  PressureIntg penalty_integrand (this);

FEM::smoothField (*p, penalty_integrand);
}
void NsPenalty1::solveAtThisTimeStep ()
{
fillEssBC ();
// initialisiert die nichtlineare Lösung mit der Lösung des
// vorherigen Zeitschrittes:
dof->field2vec (*u_prev, nonlin_solution);
// Aktualisierung der Randbedingungen:
dof->insertEssBC (nonlin_solution);

#ifdef DP_DEBUG
nonlin_solution.print(s_o,"nonlin_solution before
                        nonlinear iteration");
#endif
if (!nlsolver->solve ())
  errorFP("NsPenalty1::SolveAtThisTimeStep",
          "The nlsolver.solve call: divergence of solver \"%s\"",
          nlsolver_prm->method.c_str());
// In jedem Iterationsschritt der Schleife  nlsolver->solve,
// wird makeAndSolveLinearSystem aufgerufen!
dof->vec2field (nonlin_solution, *u);
                    // weist Lösung u zu
#ifdef DP_DEBUG
for (int k = 1; k <= grid->getNoSpaceDim (); k++)
  u()(k).values().print(s_o,oform("u(%d) after
                           nonlinear it.",k));
#endif
}
void NsPenalty1:: makeAndSolveLinearSystem ()
{
dof->vec2field (nonlin_solution, *u);
if (nlsolver->getCurrentState().method

                               == NEWTON_RAPHSON)
  dof->fillEssBC2zero();
else
  dof->unfillEssBC2zero();
makeSystem (*dof, *lineq);
lineq->solve();
#ifdef DP_DEBUG
linear_solution.print(s_o,"linear_solution");
#endif
}
void NsPenalty1:: calcElmMatVec
(int e, ElmMatVec& elmat, FiniteElement& fe)
{
// itg_rules is inherited from base class FEM
itg_rules.setRelativeOrder (0);
fe.refill (e, itg_rules);
integrands_tp = ORDINARY_TERMS;
numItgOverElm (elmat, fe);
                // Integration der lambda-freien Terme
itg_rules.setRelativeOrder (-1);
fe.refill (e, itg_rules);
integrands_tp = LAMBDA_TERMS;
numItgOverElm (elmat, fe);
               // "reduzierte" Integration der lambda-Terme
}
void NsPenalty1:: integrands
(ElmMatVec& elmat, const FiniteElement& fe)
{
if (integrands_tp == LAMBDA_TERMS)
  integrandsReduced (elmat, fe);
else if (integrands_tp == ORDINARY_TERMS)
  integrandsNonReduced (elmat, fe);
}
void NsPenalty1::integrandsNonReduced
(ElmMatVec& elmat, const FiniteElement& fe)
{

int i,j;    // Indizes der Basisfunktionen N(i)
int k,r,s;  // 1,..,nsd Indizes für
            // räumliche Dimensionen
int ig,jg;  // Zuordnungsindex für Elementmatrizen
            // und -vektoren,
            // basierend auf i,j,r,s
int e, jx;
real nabla2,h1,h2,h3,h4;
const int nsd = fe.getNoSpaceDim();
const int nbf = fe.getNoBasisFunc();
const real detJxW = fe.detJxW();
for (k = 1; k <= nsd; k++) {
  gradu_pt(k).redim (nsd);  gradup_pt(k).redim (nsd);
  u_pt(k) = u()(k).valueFEM (fe);
  u()(k).derivativeFEM (gradu_pt(k), fe);
}
real dirchl,con,dij,cij,sum;
e = fe.getElmNo();
if (nlsolver->getCurrentState().method
                               == NEWTON_RAPHSON)
  {
    for (i = 1; i <= nbf; i++) {
      for (j = 1; j <= nbf; j++) {
        jx = grid->loc2glob(e,j);
        nabla2 = 0;
        for (k = 1; k <= nsd; k++)
          nabla2 += fe.dN(i,k)*fe.dN(j,k);
        dirchl = mu*nabla2;
        sum = 0;
        for (k = 1; k <= nsd; k++)
          sum += fe.dN(j,k)*u_pt(k);
        con = density*fe.N(i)*sum;
        dij = dirchl + con;

cij = density*fe.N(i)*fe.N(j);
for (r = 1; r <= nsd; r++)
  for (s = 1; s <= nsd; s++) {
    ig = nsd*(i-1)+r;
    jg = nsd*(j-1)+s;
    h1 = cij*gradup_pt(r)(s);
    elmat.A(ig,jg) += h1*detJxW;
    if (r == s) {
      h2 = dij;
      elmat.A(ig,jg) += h2*detJxW;
      real sli = u()(r).values()(jx);
      h3 = dij*sli;
      elmat.b(ig) += -h3*detJxW;
    }
  }
}
}
}
else if (nlsolver->getCurrentState().method
                            == SUCCESSIVE_SUBST)
{
for (i = 1; i <= nbf; i++) {
for (j = 1; j <= nbf; j++) {
jx = grid->loc2glob(e,j);
nabla2 = 0;
for (k = 1; k <= nsd; k++)
  nabla2 += fe.dN(i,k)*fe.dN(j,k);
dirchl = mu*nabla2;
sum = 0;
for (k = 1; k <= nsd; k++)
  sum += fe.dN(j,k)*u_pt(k);
con = density*fe.N(i)*sum;

        dij = dirchl + con;
        cij = density*fe.N(i)*fe.N(j);

        for (r = 1; r <= nsd; r++)
          for (s = 1; s <= nsd; s++) {
            ig = nsd*(i-1)+r;
            jg = nsd*(j-1)+s;
            if (r == s) {
              h1 = dij;
              elmat.A(ig,jg) += h1*detJxW;
              h2 = 0;
              elmat.b(ig) += 0;
            }
          }
      }
    }
  }
else
  fatalerrorFP("NSPenalty1::integrands",
               "current nonlinear method=%d, illegal value",
               nlsolver->getCurrentState().method);
}

void NsPenalty1:: integrandsReduced
(ElmMatVec& elmat, const FiniteElement& fe)
{
int i,j;      // Indizes der Basisfunktionen N(i)
int k,r,s;    // 1,..,nsd Indizes für
              // räumliche Dimensionen
int ig,jg;    // Zuordnungsindex für Elementmatrizen
              // und -vektoren,
              // basierend auf i,j,r,s
real h1,h2;
const int nsd = fe.getNoSpaceDim();
const int nbf = fe.getNoBasisFunc();

const real detJxW = fe.detJxW();

for (k = 1; k <= nsd; k++) {
  gradu_pt(k).redim (nsd);   gradup_pt(k).redim (nsd);
  u_pt(k) = u()(k).valueFEM (fe);
  u()(k).derivativeFEM (gradu_pt(k), fe);
}
real sum, div;
if (nlsolver->getCurrentState().method
                               == NEWTON_RAPHSON)
  {
    for (i = 1; i <= nbf; i++)
      for (j = 1; j <= nbf; j++)
        for (r = 1; r <= nsd; r++)
          for (s = 1; s <= nsd; s++) {
            ig = nsd*(i-1)+r;
            jg = nsd*(j-1)+s;
            h1 = lambda*fe.dN(i,r)*fe.dN(j,s);
            elmat.A(ig,jg) += h1*detJxW;
          }
    sum = 0;
    for (k = 1; k <= nsd; k++)
      sum += gradu_pt(k)(k);
    div = lambda*sum;
    for (i = 1; i <= nbf; i++)
      for (r = 1; r <= nsd; r++) {
        ig = nsd*(i-1)+r;
        h2 = div*fe.dN(i,r);
        elmat.b(ig) += - h2*detJxW;
      }
  }
else if (nlsolver->getCurrentState().method
                                    == SUCCESSIVE_SUBST)
 {

     for (i = 1; i <= nbf; i++)
       for (j = 1; j <= nbf; j++)
         for (r = 1; r <= nsd; r++)
           for (s = 1; s <= nsd; s++) {
             ig = nsd*(i-1)+r;
             jg = nsd*(j-1)+s;
              h1 = lambda*fe.dN(i,r)*fe.dN(j,s);
              elmat.A(ig,jg) += h1*detJxW;
              h2 =0;
              elmat.b(ig) += 0;
           }
    }
  else
    fatalerrorFP("NSPenalty1::integrands",
                 "current nonlinear method=%d, illegal value",
                 nlsolver->getCurrentState().method);
}

B.3 main.cpp
#include <NsPenalty1.h>
main (int argc, const char* argv[])
{
  initDIFFPACK (argc, argv);
  global_menu.init ("Navier-Stokes", "NsPenalty1");
  NsPenalty1 simulator;
  global_menu.multipleLoop (simulator);
}

B.4 Generierung der Super-Elemente-Gitter

Zu Generierung eines Super-Elemente-Gitters benötigt man ein .geom-File zu Festlegung der Koordinaten jedes einzelnen Super-Elementes und ein .part-File zur Eingabe von Elementtyp, Feinheit der Unterteilung (divisions) und eventuellem Dehnungsfaktor (grading). Als Beispiel werden im folgenden die entsprechenden Files für den Fall mit 75^{∘} Öffnungswinkel vorgestellt.

Die Eingaben in supels75.geom haben die folgende Gestalt:

>no_of_dimensions=2; >subdomains=1; >no_of_supels=22;
>no_of_ind=5; >name links rechts oben unten klappe;
!Element 1
>SupEl; >subdomain_no=1; >elementtype = ElmB8n2D;
>boundary = [1(3)][3(2)];
>nodes = [1(0 0.0102)]+[2(0.006 0.0102)]+[3(0 0.02)]
          +[4(0.006 0.0165)]+[6(0.0035 0.018)];
>sides = ;
!Element 2
>SupEl; >subdomain_no=1; >elementtype = ElmB8n2D;
>boundary = [1(3)];
>nodes = [1(0 0.00281)]+[2(0.006 0.00281)]+[3(0 0.0102)]
          +[4(0.006 0.0102)];
>sides = [2(1 4)];
!Element 3
>SupEl; >subdomain_no=1; >elementtype = ElmB8n2D;
>boundary = [5(1)][1(3)];
>nodes = [1(0 0.00194)]+[2(0.006 0.00194)]+[3(0 0.00281)]
          +[4(0.006 0.00281)];
>sides = [2(2 4)];

usw. bis zu Super-Element Nr. 22. Das zugehörige supels75.part-File beginnt wie folgt:

>nsd=2; >no_of_supels=22;
!Element 1
>SupEl; >nsd=2;
>elementtype = ElmB4n2D; >divisions = [10,10]; >grading = [1,1];

!Element 2
>SupEl; >nsd=2;
>elementtype = ElmB4n2D; >divisions = [10,8]; >grading = [1,1];
!Element 3
>SupEl; >nsd=2;
>elementtype = ElmB4n2D; >divisions = [10,2]; >grading = [1,1];

Der Befehl

makegrid +iscl -c supels75 -m PreproSupElSet
               -g FILE=supels75.geom -p FILE=supels75.part -d

konstruiert dann das gesuchte Gitter. (supels75.grid)

B.5 Input-File supels75.i für Super-Elemente-Gitter

set gridfile = Verify/supels75.grid
set redefine boundary indicators = nb=4 names=inlet u0 v0 w0
                                1=(1) 2=(3 5) 3=(3 4 5) 4=(3 4 5)
set add boundary nodes = NONE
set viscosity = 4.0e-3
set penalty parameter = 800
set density = 1000
set inlet profile = 1
set characteristic inlet velocity = 20.0e-2
! 30 würde Re=1500 liefern....
set time integration parameters = dt=0
!
sub NonLinEqSolver_prm
set nonlinear iteration method = SuccessiveSubst
set max nonlinear iterations = 100
set max estimated nonlinear error = 3.0e-3
set nonlinear relaxation prm = 0.7
ok
sub LinEqAdmFE

sub Matrix_prm
set matrix type = MatBand
set symmetric storage = false
ok
sub LinEqSolver_prm
set basic method = GaussElim
ok
ok
sub SaveSimRes
set field storage format = BINARY
set grid storage format = BINARY
ok
ok

B.6 Input-File king.i für BOX WITH BELL-Geometrie

set gridfile = Verify/supels75.grid
set redefine boundary indicators = nb=4 names=inlet u0 v0 w0
                                  1=(1) 2=(3 5) 3=(3 4 5) 4=(3 4 5)
set viscosity = 4.0e-3
set penalty parameter = 800
set density = 1000
set inlet profile = 1
set characteristic inlet velocity = 20.0e-2 ! 30 gives Re=1500
set time integration parameters = dt=0
!
sub NonLinEqSolver_prm
set nonlinear iteration method = SuccessiveSubst
set max nonlinear iterations = 100
set max estimated nonlinear error = 3.0e-3
set nonlinear relaxation prm = 0.7
ok
sub LinEqAdmFE
sub Matrix_prm
set matrix type = MatBand
set symmetric storage = false
ok
sub LinEqSolver_prm
set basic method = GaussElim
ok

ok
sub SaveSimRes
set field storage format = BINARY
set grid storage format = BINARY
ok
ok

Literatur

[1] J. Aagaard, C.N. Hansen, J. Tingleff, and I. Rygg. Seven-and-a-half Years Clinical Experience with the Carbomedics Prosthetic Heart Valve. J. Heart Valve Dis., 4:628–633, 1995.

[2] B.J. Bellhouse and F.H. Bellhouse. Fluid mechanics of model normal and stenosed artic valves. Circ Res, 25:693–704, 1969.

[3] B.J. Bellhouse and F.H. Bellhouse. Fluid mechanics of model mitral valve and left ventrikel. Cardiovasc Res, 6:199–210, 1972.

[4] B.J. Bellhouse and L. Talbot. J. Fluid Mech., 35:721, 1969.

[5] A. M. Bruaset. A survey of preconditioned iterative methods. Longman Scientific & Technical, 1. edition, 1995.

[6] Sulzer Carbomedics. Infomaterial, 1998.

[7] A.J. Chorin and J.E. Marsden. A Mathematical Introduction to Fluid Mechanics. Springer Verlag, 1979.

[8] L. Collatz. Differentialgleichungen. Teubner Verlag, 6. edition, 1981.

[9] C.W. Cryer. Numerik Partieller Differentialgleichungen 1. WWU Münster, 1998.

[10] C.W. Cryer. Numerik Partieller Differentialgleichungen 2. WWU Münster, 1998.

[11] C.W. Cryer and P.P. Lunkenheimer. Mathematical Modelling of the Cardiovascular System. WWU Münster, 1995.

[12] Czhiak, Langer, and Ziegler. Biologie. Springer Verlag, 5. edition.

[13] P. Deuring and W. von Wahl. Strong Solutions of the Navier-Stokes System in Lipschitz Bounded Domains. Math. Nachr., 171:111–148, 1995.

[14] O. Forster. Analysis 2. Vieweg Studium, 5. edition, 1993.

[15] Y.C. Fung. Biomechanics, Circulation. Springer Verlag, 2. edition.

[16] A. Greenbaum. Iterative Methods for Solving Linear Systems. SIAM, 1. edition, 1997.

[17] Ch. Großmann and H.-G. Roos. Numerik partieller Differentialgleichungen. Teubner Verlag, 1. edition, 1992.

[18] Y. Henderson and F. Johnson. Two models of closure of heart valves. Heart, 4:69–82, 1912.

[19] D. Horstkotte and F. Loogen. Update in Heart Valve Replacement. Proceedings of the Second European Symposium on the St.Jude Medical Heart Valve, 1986.

[20] J. Jost. Partielle Differentialgleichungen. Springer Verlag, 1. edition, 1998.

[21] R. Juchems. Herz- und Kreislaufkrankheiten. Wiss. Buchges. Darmstadt.

[22] T. Kato. Strong L^{p}-Solutions of the Navier-Stokes Equation in ℝ^{m}, with Applications to Weak Solutions. Math. Z., 192:135–148, 1986.

[23] Mary J. King. Computational and experimental studies of flow through a bileaflet mechanical heart valve. PhD thesis, University of Leeds, U.K., 1994.

[24] H. Kozono and T. Ogawa. Some L^{p} Estimate for the Exterior Stokes Flow and An Application to the Non-Stationary Navier-Stokes Equations. Indiana Univ. Math. J., 41:789–808, 1992.

[25] W. Krahwinkel, S. Moltzahn, and M. Zeydabadinejad. Echokardiographie der künstlichen Herzklappen. Thieme Verlag, 1995.

[26] H. P. Langtangen. Computational partial differential equations: Numerical methods and Diffpack Programming. Springer Verlag, 1. edition, 1999.

[27] H. P. Langtangen. Finite element prepocessors in Diffpack. Numerical Objects Series, Numerical Objects A.S., 1999.

[28] C.S.F. Lee and L. Talbot. A fluid mechanical study on the closure of heart valves. J Fluid Mech, 91:41–63, 1979.

[29] H. Leonhardt. Taschenatlas der Anatomie. Thieme Verlag, 5. edition.

[30] C.S. Peskin. Numerical Analysis of Blood Flow in the Heart. J. Computational Physics, 25:220–252, 1977.

[31] W. Richter. Partielle Differentialgleichungen. Spektrum, 1. edition, 1995.

[32] H. R. Schwarz. Methode der finiten Elemente. Teubner Verlag, 1. edition, 1980.

[33] J. S. Schwegler. Der Mensch – Anatomie und Physiologie. Thieme Verlag, 2. edition.

[34] E. M. Stein. Singular Integrals and Differentiability of Functions. Princeton University Press, 1. edition, 1970.

[35] R. Temam. Navier-Stokes Equations. North Holland Publishing, 1. edition, 1977.

[36] R. Temam. Navier-Stokes Equations and Nonlinear Functional Analysis. SIAM, 1. edition, 1983.

[37] Prof. W. Walter. Gewöhnliche Differentialgleichungen. Springer Verlag, 5. edition, 1991.

[38] K. Wieghardt. Theoretische Strömungslehre. Teubner Verlag, 2. edition, 1974.

[39] M. Wiegner. The Navier-Stokes Equations - a Neverending Challenge? Jber. d. Dt. Math.-Verein., 101:1–25, 1999.