Der erfolgreiche Erstflug einer amerikanischen Raumfähre samt weicher Landung auf einem gewöhnlichen Rollfeld war ein weltweit verfolgtes Ereignis. Dagegen wurde eine spätere Meldung der Luft- und Raumfahrtbehörde NASA kaum mehr beachtet: Die am 12. April 1981 gestartete Columbia war nur knapp einer Katastrophe entgangen.

Während der ersten Phase des Wiedereintritts in die Erdatmosphäre mußten die Astronauten die Nase des Space Shuttle um 40 Grad gegen die Flugrichtung aufrichten, indem sie eine Steuerklappe in den Fahrtwind stellten. Abschätzungen hatten zuvor ergeben, daß eine Auslenkung um 7 Grad aus der Ruhelage den beabsichtigten Effekt haben sollte. Das erwies sich jedoch in der Realität als völlig unzureichend. Nur der erhebliche Sicherheitszuschlag, den die Konstrukteure vorgesehen hatten, rettete der Besatzung das Leben: Erst durch Ausfahren bis auf 16 Grad brachten die Piloten das Raumfahrzeug in die richtige Position; bei 21 Grad war der Anschlag.

Die Ursachenforschung ergab, daß die Konstrukteure einen Effekt nicht richtig berücksichtigt hatten, der sich erst bei den extrem hohen Geschwindigkeiten des Wiedereintritts einstellt: Die Luft wird durch einen mit 24facher Schallgeschwindigkeit eindringenden Flugkörper so stark komprimiert und erhitzt, daß die Sauerstoff- und Stickstoffmoleküle in einzelne Atome aufbrechen. Dadurch änderten sich bei der Rückkehr der Columbia die thermodynamischen Eigenschaften der Luft und damit wiederum die Wirksamkeit der Steuerklappe.

Mittlerweile ist dieser Effekt zwar im Prinzip kalkulierbar; es verbleiben jedoch Unsicherheiten. Wesentliche Größen wie etwa die Reaktionsraten der maßgeblichen chemischen Prozesse sind im Labor nur schwer zu bestimmen; Extrapolationen aus den Verhältnissen bei weniger extremen Temperaturen sind mit großen Unsicherheiten behaftet. Diese schlagen durch bis auf die Endergebnisse einer Computersimulation, so daß man sie durch geeignete Experimente absichern muß.

Die Deutsche Forschungsanstalt für Luft- und Raumfahrt (DLR) hat in den vergangenen knapp zehn Jahren im Rahmen der europäischen Raumfahrtprogramme die erforderliche Symbiose zwischen Berechnung und Experiment besonders intensiv gepflegt: einerseits mit einer speziell gebauten Versuchsanlage, dem Hochenthalpie-Stoßkanal im DLR-Forschungszentrum Göttingen (siehe Kasten auf der letzten Doppelseite) , andererseits durch die Entwicklung spezialisierter numerischer Verfahren. In diesem Artikel geht es vor allem um die numerische Seite des Problems.

Obgleich die weitere Entwicklung der europäischen Raumfähre Hermes mittlerweile nicht mehr betrieben wird, sind unsere Forschungen keineswegs gegenstandslos geworden. Insbesondere wären das von europäischer Seite geplante Crew Transfer Vehicle (CTV) für den Pendelverkehr zur Raumstation Alpha oder die von der NASA vorgeschlagenen Aeroassisted Space Transfer Vehicles (ASTV; Bild 2) ohne numerische Strömungssimulation und entsprechende Windkanalexperimente nicht zu bauen. Ein ASTV ist ein wiederverwendbares Raumfluggerät, das beim Eintritt in eine Planetenatmosphäre oder beim Übergang von einer hohen in eine niedrigere Umlaufbahn seine Geschwindigkeit nicht – wie bisher üblich – durch Raketenantriebe, sondern allein durch das Ausnutzen aerodynamischer Kräfte vermindert. Das erspart erhebliche Mengen Treibstoff; entsprechend läßt sich die Nutzlast erhöhen. Mit einem ASTV wären Reparatur oder Ausbau teurer Kommunikations- und Forschungssatelliten nach Einschätzung der NASA weit kostengünstiger als Bau und Aussetzen eines neuen Systems. Ein bemannter Flug zum Mars gerät erst durch das ASTV-Prinzip mit seinem verminderten Treibstoffbedarf in den Bereich des Möglichen.


Verdünnte und Kontinuumsströmung

Wenn ein Flugkörper wieder in die Erdatmosphäre eintritt, gerät er zunächst in sehr dünne Luft. Indem die Gasmoleküle auf seine Hülle prasseln, entziehen sie ihm Bewegungsenergie und wandeln sie in Wärme um. Aber sie treten nicht nennenswert miteinander in Wechselwirkung, so daß man von makroskopischen Größen wie Druck und Temperatur der Luft nicht sprechen kann.

Das beeinflußt die mathematische Modellierung. Im Prinzip müßte man Bewegungsgleichungen für jedes einzelne Luftmolekül lösen. Weil das für zunehmende Teilchenzahlen in dichter werdender Atmosphäre nicht praktikabel ist, verwendet man eine statistische Mittelung, was auf die sogenannte Boltzmann-Gleichung hinausläuft (vergleiche Spektrum der Wissenschaft, Februar 1995, Seite 21); der österreichische Mathematiker und Physiker Ludwig Boltzmann (1844 bis 1906) hat die statistische Thermodynamik begründet. Diese – sehr komplizierte – Gleichung ist eine mathematisch korrekte Beschreibung der physikalischen Verhältnisse; aber ihre numerische Lösung ist ebenfalls nicht sonderlich praktikabel. Statt dessen läßt man in der Computersimulation stellvertretend einige zehntausend statt hundert Billionen Moleküle fliegen; die Statistik sorgt dann für das korrekte Ergebnis – vorausgesetzt, überall im Rechenraum sind genug Stellvertreter-Moleküle für eine statistische Mittelung vorhanden. Dieses Verfahren heißt direkte Monte-Carlo-Simulation (Spektrum der Wissenschaft, April 1994, Seite 64).

Der Flug in diesen großen Höhen ist zwar schwierig zu modellieren und zu berechnen, bereitet aber dem Ingenieur nicht die größte Sorge. Kritisch ist vielmehr jener Bereich der Flugbahn, in dem sich der Flugkörper am meisten aufheizt. Auch eine theoretische Bahn für ihn zu berechnen, die am Boden mit der Geschwindigkeit null endet, ist kein großes Problem. Schwer wird es erst durch das, was im Formalismus eine Nebenbedingung heißt: Die Temperatur darf einen gewissen Wert nicht überschreiten, da-mit Raumfahrzeug und Insassen keinen Schaden nehmen. Die maximale Wärmelast tritt bei der in Bild 3 angegebenen Bahn in einer Höhe von etwa 70 Kilometern auf. Dort ist die Luft schon etwa ein zehntausendstel so dicht wie in Meereshöhe; sie läßt sich für unsere Zwecke als Kontinuum auffassen, so daß Dichte, Druck und Temperatur in jedem Punkt wohldefinierte Größen sind.

Es geht also darum, detailliert zu berechnen, wie die Luft den Flugkörper umströmt. Dazu muß man wissen, welche physikalischen und chemischen Eigenschaften sie in der Umgebung des Körpers annimmt. Daraus ergibt sich schließlich, wie stark dieser abgebremst und erwärmt wird.

Zweckmäßig wählt man ein Koordinatensystem, das mit dem Flugkörper mitfliegt; man fragt also nach dem Strömungsbild, das sich einstellt, wenn ein ruhender Flugkörper von vorn angeblasen wird. Wie bei anderen strömungsmechanischen Problemen üblich, beschreibt man das physikalische System durch die Funktionen Dichte, Geschwindigkeit, Temperatur und Druck, die an jedem Punkt des Rechengebietes einen gewissen Wert annehmen. Sind diese Werte zu einem gewissen Zeitpunkt sämtlich bekannt, so gibt ein System partieller Differentialgleichungen, das der Navier-Stokes-Gleichungen, Auskunft darüber, wie sich diese Werte in dem nächsten kleinen Zeitintervall verändern werden. Somit besteht die Aufgabe darin, diese Gleichungen, die der französische Ingenieur Claude Navier (1785 bis 1838) und der englische Physiker und Mathematiker Sir George Gabriel Stokes (1819 bis 1903) aufgestellt haben, zu lösen – mit den für diesen Fall typischen Erschwernissen.

Erstens schiebt der Flugkörper, wie jeder mit Überschallgeschwindigkeit bewegte Körper, eine Stoßwelle vor sich her (wenn diese Welle unser Ohr trifft, hören wir den bekannten Überschallknall). An der Stoßfront ändern sich die thermodynamischen Variablen sprunghaft: Die Temperatur steigt von einem sehr niedrigen Wert schlagartig auf mehrere tausend Kelvin, die Anströmungsgeschwindigkeit fällt plötzlich von 7 Kilometer pro Sekunde auf ungefähr ein Hundertstel dieses Wertes, und für die anderen Größen gilt ähnliches (Bild 2).

Das ist schon für die Theorie ein Problem, denn die bloße Gültigkeit einer Differentialgleichung wie der von Navier und Stokes impliziert, daß ihre Variablen bezüglich des Raumes stetig sind, also gerade keine Sprünge machen. Mehr noch – die meisten numerischen Verfahren sind auf stetige Lösungen ausgelegt und versagen, wenn eine unstetige zu finden ist. Es sind also spezielle Anpassungen vorzunehmen.

Zweitens wird durch Reibungskräfte eine Schicht von Luftmolekülen an der Oberfläche festgehalten; aber die Wirkung dieser Randbedingung reicht bei Hyperschallströmungen nicht weit in die umgebende Luft hinaus. Die Verhältnisse ändern sich also im allgemeinen quer zur Oberfläche viel dramatischer als entlang derselben, was bei einigen numerischen Verfahren ebenfalls Schwierigkeiten bereitet.

Physikalisch-chemische Vorgänge

Drittens werden durch den starken Temperaturanstieg hinter dem Verdichtungsstoß chemische Reaktionen ausgelöst. Die Moleküle geraten in heftige innere Vibrationen und brechen schließlich auseinander. Da diese Dissoziation Energie verzehrt und deswegen die Temperatur mindert, ist sie im Prinzip sehr erwünscht. Bei sehr hohen Temperaturen werden die Atome auch noch ionisiert.

Bis eine solche Zustandsänderung vollzogen ist, vergeht indes eine gewisse Zeit, weil nur durch Zusammenstöße unter den Molekülen Bewegungsenergie sich in innere Energie verwandeln kann. Je größer die Höhe über Grund, desto dünner ist die Luft und – wegen der geringeren Kollisionsrate – desto länger diese Relaxationszeit. In großen Höhen kann es sein, daß die Reaktion kaum begonnen hat, wenn die Luft schon längst am Flugkörper vorbeigeströmt ist.

Das Verhältnis zwischen der Zeitspanne, in der ein Partikel einen Weg von der Gesamtlänge des Fluggeräts zurücklegt, und der Relaxationszeit wird als Damköhler-Zahl bezeichnet. (Der deutsche Chemiker Gerhard Damköhler, 1908 bis 1944, war der Begründer der chemischen Reaktionstechnik.) Im Grenzfall einer sehr großen Damköhler-Zahl ist die Relaxationszeit so kurz, daß sich praktisch sofort ein Gleichgewichtszustand einstellt. Bei sehr kleinen Damköhler-Zahlen findet während des interessierenden Zeitraums praktisch keine chemische Reaktion statt. Man spricht dann von einem eingefrorenen Zustand und entsprechend von eingefrorener Strömung (frozen flow) – bei Temperaturen von mehreren tausend Grad eine einigermaßen paradoxe Bezeichnung.

Für den kritischen Teil der Wiedereintrittsbahn ist allerdings die Damköhler-Zahl weder sehr klein noch sehr groß; der chemische Zustand der Luft ändert sich also, während sie am Flugkörper vorbeiströmt. Die Beschreibung dieses Zustandes erfordert somit außer den üblichen Variablen noch zahlenmäßige Angaben darüber, wie die Luft an jedem Punkt des Rechengebiets chemisch zusammengesetzt ist: Außer molekularem Stickstoff und Sauerstoff kommen beide Elemente auch atomar vor, dazu kompliziertere Verbindungen wie Stickoxide. Die Navier-Stokes-Gleichungen sind entsprechend zu erweitern um Angaben, wie sich der Anteil jeder Stoffsorte (Spezies) durch Strömung und chemische Reaktionen mit der Zeit verändert.

Die Luft kann sich des weiteren nicht nur chemisch, sondern auch thermisch im Nichtgleichgewicht befinden. Die Moleküle haben dann nicht ausreichend Zeit, ihre Vibrations- und sonstigen Bewegungsenergien durch Kollisionen so untereinander auszutauschen, daß sich eine einheitliche Temperatur einstellt. In diesem Falle muß man die Vibrationsenergie der Moleküle jeder Spezies als weitere Variable einführen.

Insgesamt ist die numerische Behandlung von Strömungen im chemischen und thermischen Nichtgleichgewicht wesentlich mühsamer als die von eingefrorenen oder Gleichgewichtsströmungen. Dieser Aufwand ist jedoch unvermeidlich, weil die vereinfachte Modellierung – zum Beispiel für die Wirksamkeit der Steuerklappe – nur unbefriedigende Ergebnisse liefert.

Schließlich schlägt die Strömung in der wandnahen Schicht, der Grenzschicht des Raumflugzeugs, in späteren Phasen des Wiedereintritts von laminar zu turbulent um. Statt auf geordneten, glatten und damit vorhersagbaren Bahnen bewegen sich die Teilchen dann äußerst ungeordnet, weil Wirbel aller Größenordnungen auftreten.


Navier-Stokes-Gleichungen

Die Umsetzung der physikalischen Gesetze in den mathematischen Formalismus folgt einem einfachen Rezept: Es gelten die Erhaltungssätze für Masse, Impuls und Energie. Man denke sich ein kleines quaderförmiges Testvolumen irgendwo im Rechengebiet. Die darin enthaltene Gasmasse kann sich nur dadurch verändern, daß Masse in das Testvolumen hinein- oder aus ihm herausströmt. Für Impuls und Energie gilt der gleiche Gedankengang. Die Navier-Stokes-Gleichungen ergeben sich, indem man Bilanzen dieser physikalischen Größen aufstellt und dann das Testvolumen durch einen Grenzwertprozeß zu einem Punkt zusammenschrumpfen läßt.

Genauer besehen, kann die jeweilige Größe durch den Wind (konvektiv) transportiert werden oder sich durch jene ungerichtete Bewegung (diffusiv) ausbreiten, die für die Wärmeleitung typisch ist. Die Bilanzänderung durch Konvektion ist proportional der räumlichen Änderung der betreffenden Größe: Die Temperatur nimmt an einem bestimmten Ort durch Konvektion nur zu, wenn es dort, wo der Wind herkommt, heißer ist. Mathematisch wird das durch eine Ableitung (den Differentialquotienten) nach dem Ort ausgedrückt. Bilanzänderung durch Diffusion hingegen ist der zweiten Ableitung nach dem Ort proportional: An einem linearen Temperaturgefälle ändert die Wärmediffusion nichts, wohl aber flacht sie lokale Spitzen ab und füllt lokale Senken auf.

Schließlich gibt es Bilanzänderungen, die nicht von den Verhältnissen in der Umgebung des geschrumpften Testvolumens abhängen, sondern nur von den Werten in dem Punkt selbst. Dazu gehören die Änderungen der Massenanteile der einzelnen Spezies durch chemische Reaktionen.

Insgesamt ergibt sich, da auch noch alle drei Raumrichtungen zu berücksichtigen sind, ein recht kompliziertes System von Gleichungen, die als nichtlineare partielle Differentialgleichungen zweiter Ordnung klassifiziert werden. Für die Zwecke dieses Artikels genügt es jedoch, der Einfachheit halber nur eine Raumdimension und eine abhängige Variable (beispielsweise die Temperatur) zu betrachten. Es ergibt sich eine schematische Gleichung (Bild 4 oben), an der man bereits die wesentlichen Probleme zeigen kann.

Die drei physikalischen Prozesse Konvektion, Diffusion und Reaktion haben sehr verschiedene Auswirkungen. Die Konvektion transportiert den Zustand u, der in einem gewissen Punkt des Rechengebiets besteht, mit der Geschwindigkeit a(u) weiter (vergleiche die Gleichung in Bild 4). Im Raum-Zeit-Diagramm ergibt sich eine Kurve, die sogenannte Charakteristik. Aber die Geschwindigkeit a ist ihrerseits vom Zustand u abhängig, also im allgemeinen nicht räumlich konstant. Deswegen kann es vorkommen, daß zwei (gleichläufige) Charakteristiken sich schneiden. Im Schnittpunkt müßte also das System zwei verschiedene Zustände annehmen, was unmöglich ist. Statt dessen bildet sich eine Unstetigkeit: ein Punkt, zu dessen beiden Seiten der Systemzustand deutlich verschiedene Werte annimmt. Ein klassisches Beispiel dafür ist die erwähnte Überschall-Stoßfront.

Bei der Diffusion hingegen wird Information nicht entlang ausgezeichneter Kurven transportiert, sondern ist sofort überall im Raum spürbar. Während Konvektionseffekte tendenziell Unterschiede verschärfen, wirkt die Diffusion als großer Gleichmacher, und das um so intensiver, je größer der Koeffizient e ist.

Der Koeffizient e ist als Kehrwert einer dimensionslosen Kennzahl, der Reynolds-Zahl, zu interpretieren. Kleine Werte der Reynolds-Zahl (große von e) korrespondieren mit zähen Strömungen, in denen die Reibung überwiegt; diese Strömungen sind vorwiegend laminar. Für größere Reynolds-Zahlen (kleinere Werte von e) überwiegen mehr und mehr die Trägheitskräfte. Auch in den Grenzschichten kann dann der Dämpfungseffekt der Reibung nicht mehr das Auftreten von Turbulenz verhindern. Deren Struktur reicht jedoch bis in so geringe Dimensionen, daß sie auch mit heutiger Rechnerkapazität nur in unrealistisch einfachen Fällen direkt berechenbar ist. Man muß sich weiterhin damit behelfen, die Navier-Stokes-Gleichungen – durch die sogenannte Favre-Mittelung – zu modifizieren und die makroskopischen Auswirkungen der Turbulenz, vor allem die zusätzlich erzeugte Zähigkeit, durch ein geeignetes Modell näherungsweise zu erfassen.

Die Prozesse der dritten Art, die wir unter "Reaktion" zusammengefaßt haben, sind in der Theorie nicht so problematisch wie die beiden anderen. Allerdings laufen die chemischen Reaktionen mit höchst unterschiedlichen Geschwindigkeiten ab. Differentialgleichungen mit der zugehörigen Eigenschaft nennt man steif. In solchen Systemen dominieren die kurzzeitigen, starken Änderungen alle anderen Phänomene, woran das Rechenverfahren angepaßt sein muß.


Numerische Methoden

Die genannten physikalischen Prozesse sind von so unterschiedlicher Art, daß man sie zweckmäßigerweise mit völlig verschiedenen numerischen Verfahren angeht. Man kann sogar diese Vorgänge, die alle gleichzeitig stattfinden, nacheinander behandeln – so, als würde auf einen vorliegenden Zustand für eine gewisse kurze Zeitspanne Dt nur die Konvektion einwirken, dann für dieselbe Zeitspanne nur die Diffusion und dann nur die Reaktion (Bild 4 unten). Das Verfahren heißt operator splitting, weil man den Operator (das mathematische Gebilde), der zu einem Zustand dessen zeitliche Änderung bestimmt, in drei Teile aufspaltet.

Dabei führt man allerdings einen Fehler in die Berechnung ein; er ist jedoch qualitativ nicht schlimmer als der, den man ohnehin machen muß: Die korrekte Gleichung setzt die zeitliche Änderung im Grenzwert unendlich kleiner Zeiträume (¦u/¦t) mit der räumlichen Änderung im Grenzwert unendlich kleiner Entfernungen (¦u/¦x) in Beziehung. Für die numerische Berechnung muß man das Problem diskretisieren, das heißt zumindest bezüglich der Zeit die genannten Differentialquotienten – das sind Grenzwerte von Differenzenquotienten – durch geeignete Differenzenquotienten ersetzen. Operator splitting verursacht einen Fehler in der Größenordnung des zeitlichen Diskretisierungsfehlers.

Was den Raum betrifft, so steht man vor demselben prinzipiellen Problem wie bezüglich der Zeit: Es ist eine Gleichung für Funktionen zu lösen, die in unendlich vielen Punkten einen Wert haben und bei denen es auf Unterschiede dieser Werte in beliebig eng benachbarten Punkten ankommt. Man kann aber nur mit endlich vielen Unbekannten rechnen; wenn das Funktionswerte in gewissen Punkten sein sollen, können diese nicht beliebig eng benachbart sein. In der Strömungsmechanik sind im wesentlichen drei Näherungsverfahren gebräuchlich (Bild 5).

Das theoretisch einfachste ist das der finiten Differenzen (FD), das auch für die zeitliche Diskretisierung angewandt wird: Man überzieht das Rechengebiet mit einem Rechtecksgitter und setzt Differenzenquotienten zwischen benachbarten Gitterpunkten an die Stelle der Ableitungen. (Es ist möglich und häufig sinnvoll, das Gitter in Anpassung an die Geometrie des Problems geeignet zu verzerren.) Durch diese Ersetzung macht man einen Fehler beim Aufstellen der Gleichungen, die man hinterher löst. Je feiner das Gitter, desto kleiner ist dieser Fehler. Nur gibt es keinen einfachen Zusammenhang zwischen dem Fehler der Gleichung und dem der Lösung dieser Gleichung; das ist – außer der Beschränkung auf die regelmäßige Gitterform – die wesentliche Schwäche der FD-Verfahren.

Dagegen diskretisieren die Finite-Elemente-(FE)-Verfahren nicht das Problem, sondern die (unbekannte) Lösung: Man approximiert die Lösung der Differentialgleichung durch eine Summe einfacherer Funktionen (meist Polynome) und fordert, daß diese, in die Gleichung eingesetzt, einer Lösung so nahe kommt wie überhaupt möglich. Daraus gewinnt man gewöhnliche Gleichungen für die Parameter dieser Funktionen. FE-Methoden sind auf beliebigen Rechennetzen einsetzbar; in der Praxis bevorzugt man Gitter, die aus Dreiecken und/oder Vierecken bestehen.

Die Finite-Volumen-(FV)-Methoden schließlich greifen den Gedanken wieder auf, der bei der Herleitung der Navier-Stokes-Gleichungen eine Rolle gespielt hat: Man teilt das Rechengebiet in kleine Teilvolumina. Die Unbekannten des Systems sind (unter anderem) Masse, Impuls und Energie, die in jedem Teilvolumen enthalten sind, und die Gleichungen unter ihnen sind Bilanzbedingungen. Es muß zum Beispiel genau so viel Masse durch die linke Wand in ein Teilvolumen einströmen, wie durch die rechte des linken Nachbarvolumens ausströmt. Dadurch ergibt sich nahezu automatisch, daß die Gesamtmasse des diskretisierten Systems erhalten bleibt; gleiches gilt für die anderen Erhaltungsgrößen wie Impuls und Energie. Das genäherte System reproduziert also eine wesentliche physikalische Eigenschaft des echten; das gilt für die beiden anderen Verfahren nur näherungsweise.

Wie die FE-Methoden sind FV-Verfahren nicht zwingend an reguläre Rechengitter gebunden. Man hat deshalb große Freiheiten, das Netz der Punkte, an denen man die Lösung berechnet, der Geometrie des Problems anzupassen. Besonders interessant ist es, zur Erhöhung der Genauigkeit zusätzliche Punkte dort einzufügen, wo die Lösung heftig variiert, und zur Einsparung überflüssiger Arbeit Punkte wegzulassen, wo sich nichts Wesentliches abspielt (Bild 8).

Dazu muß man zumindest im nachhinein (a posteriori) den berechneten Werten eine Schätzung dafür entnehmen können, wie groß der numerische Fehler im Einzelfall ist. Bei einem zu großen Fehler ist das Gitter zu verfeinern, bei einem sehr kleinen kann man es vergröbern. Für unsere Algorithmen haben wir solche Schätzverfahren in Zusammenarbeit mit Endre Süli von der Universität Oxford und Gerald Warnecke von der Universität Magdeburg entwickelt.

Da sich, wie erwähnt, die Verhältnisse quer zur Grenzschicht viel dramatischer ändern als längs, produziert ein adaptiver Algorithmus entlang der Grenzschicht extrem langgestreckte Drei- oder Vierecke. Das ist im Prinzip physikalisch sinnvoll, wirkt sich jedoch nachteilig auf eine – unten näher erläuterte – Nachbearbeitung zur Genauigkeitssteigerung aus. Abhilfe schafft eine geeignete Modifikation der randnahen Gitterzellen.


Die Lösung der Transportgleichung

Wie setzt man nun auf der Basis eines dieser Diskretisierungsverfahren die mathematische Formulierung der drei Prozesse Konvektion, Diffusion und Reaktion in jeweils eigene numerische Schemata um? Die meisten Untersuchungen gibt es zur Konvektion; denn dieser konzeptuell einfache Vorgang bietet numerisch besondere Schwierigkeiten, weil Unstetigkeiten auftreten können.

Eine der wesentlichen Ideen besteht darin, den Stier bei den Hörnern zu packen: Wenn man die Unstetigkeiten schon nicht vermeiden kann, dann tue man gleich so, als enthalte der Zustand lauter Unstetigkeiten und sei zwischen ihnen konstant – eine Vorstellung, die besonders gut zur FV-Diskretisierung paßt. Für diesen Spezialfall gibt es nämlich eine explizite Lösung des Problems. Sie geht auf den in Göttingen wirkenden Mathematiker Bernhard Riemann (1826 bis 1866) zurück, dessen Name ansonsten eher mit Gegenständen der reinen Mathematik wie dem Riemann-Integral und Riemannschen Flächen verknüpft ist.

Riemann betrachtete ein Gas in einer langen Röhre, in deren Mitte sich eine Membran befindet. Links und rechts davon sind Druck, Dichte und Temperatur verschieden. Wenn plötzlich die Membran verschwindet, entsteht eine interessante Strömung. Obwohl Riemanns Analyse dieses Vorgangs noch fehlerhaft war, nennt man solche Systeme, die tatsächlich für Experimente eingesetzt werden, Riemannsche Stoßrohre und die Aufgabe, das Verhalten des Gases mathematisch zu beschreiben, ein Riemann-Problem. Der Hochenthalpie-Stoßkanal der DLR in Göttingen ist ein Beispiel für ein (sehr modernes) Stoßrohr.

Nach dem Entfernen der Membran bewegt sich ein Verdichtungsstoß in das Gas mit dem geringeren Druck. Zurück läuft eine Verdünnungswelle, in der sich die Strömungsgrößen stetig ändern. Eine sogenannte Kontaktunstetigkeit trennt Verdünnungswelle und Verdichtungsstoß. Lediglich Druck und Geschwindigkeit sind über die Kontaktunstetigkeit hinweg stetig; alle anderen Strömungsgrößen ändern sich sprunghaft (Bild 6).

Bei gegebenem Anfangszustand läßt sich der Zustand des Systems zu jeder späteren Zeit berechnen. Diese Idee liefert die sogenannten Riemann-Löser, die Bestandteil zahlreicher moderner FD- und FV-Verfahren geworden sind: Man approximiere den Anfangszustand des Systems durch stückweise konstante Funktionen mit Sprüngen, berechne die exakte Lösung der zugehörigen Riemann-Probleme bis zu einem Zeitpunkt, zu dem die von benachbarten Unstetigkeiten ausgehenden Wellen beziehungsweise Stöße einander noch nicht beeinflussen – bis dahin ist die Rechnung noch einfach –, approximiere das Ergebnis wieder durch stückweise konstante Funktionen mit Sprüngen, und so weiter.

Während sich Verdichtungsstöße in FD- und FV-Verfahren problemlos auf wenigen Gitterzellen darstellen lassen, muß man für Finite-Elemente-Verfahren zusätzliche Kunstgriffe anwenden, die häufig noch vom jeweiligen Problem abhängen. Darum sind FE-Verfahren in der Strömungsmechanik kompressibler Medien wenig verbreitet. Gleichwohl finden sich viele Ideen aus der Methode der finiten Elemente in FV-Verfahren wieder.

Ein anderes Prinzip hat eine nicht ganz so lange Tradition. Sie beginnt mit Arbeiten des ungarisch-amerikanischen Mathematikers John von Neumann (1903 bis 1957) im Rahmen des Manhattan-Projekts – der Entwicklung der amerikanischen Kernspaltungsbombe – in den vierziger Jahren in Los Alamos (New Mexico), basiert aber auf wesentlich älteren theoretischen Arbeiten der Schule um Richard Courant (1888 bis 1972) im Göttingen der zwanziger Jahre. Die Grundidee ist, die Unstetigkeit gewissermaßen aufzuweichen, indem man dem System das künstlich hinzufügt, was auch in der Natur häufig das Auftreten von Sprüngen verhindert: ein bißchen Zähigkeit (artificial viscosity). Man setzt also an die Stelle des echten Transportproblems ein verfälschtes; aber wenn man es geschickt macht, bleibt der dadurch eingeführte Fehler wieder in der Größenordnung des Diskretisierungsfehlers. Durch die künstliche Glättung wird aus einem Verdichtungsstoß eine starke – aber differenzierbare – Kompressionswelle; dadurch ist das Problem den üblichen Verfahren, die Glattheit voraussetzen, wieder zugänglich.

Ein drittes Konzept ist das Gegenwind-(upwind-)Verfahren, das Courant in den fünfziger Jahren gemeinsam mit Mina Rees und Eugene Isaacson am Courant-Institut der Universität New York erfand. Ihm liegt der Gedanke zugrunde, daß man bei Überschallströmungen eigentlich nur in die Richtung schauen muß, aus welcher der Wind kommt, weil aus anderen Richtungen keine Information kommen kann (während im allgemeinen die besten Verfahren diejenigen sind, die alle Richtungen genau gleich behandeln). Das setzt sich um in eine zustandsabhängige Diskretisierung. Inzwischen sind Upwind-Verfahren aus der modernen numerischen Strömungsmechanik nicht mehr wegzudenken. Sie finden sich algorithmisch in approximativen Riemann-Lösern wieder und werden im Rahmen von FD- wie auch FV-Methoden gleichermaßen genutzt.

Naiv angewandt, neigen alle klassischen Diskretisierungsverfahren dazu, an kritischen Stellen zu übersteuern wie ein überdrehter Verstärker: Winzige Störungen schaukeln sich zu Schwingungen auf, die mit dem Problem nichts zu tun haben. Es handelt sich nicht, wie beim Verstärker, um zeitliche, sondern um räumliche Schwingungen; sie können so heftig sein, daß das Verfahren völlig unsinnige Werte berechnet wie beispielsweise negative Dichten oder Drücke.

Viskosität, sei sie im Problem enthalten oder künstlich eingeführt, dämpft diese Oszillationen; aber es ist schwierig, das richtige Maß an Viskosität zu finden. Der israelische Mathematiker Ami Harten hat ein Verfahren ersonnen, das direkter an die Oszillationen anknüpft. Es arbeitet mit der Variation einer Funktion; das ist das Integral über den Betrag ihrer Ableitung, anschaulich gesprochen: die Summe aller Abstände zwischen Minima und benachbarten Maxima oder umgekehrt. Die Variation einer monoton fallenden oder höchstens konstanten Funktion wie der Dichte in Bild 6 unten ist einfach die Differenz zwischen ihren beiden Randwerten. Bei der exakten Lösung nimmt die Variation nicht zu; insbesondere können sich – wenn die Randwerte unverändert bleiben – keine neuen Maxima oder Minima bilden.

Die Gesamtvariation einer echten Lösung der Navier-Stokes-Gleichungen kann aus physikalischen Gründen nicht zu-, sondern höchstens abnehmen. Verfahren, die diese Eigenschaft in der numerischen Lösung reproduzieren, heißen TVD (total variation diminishing; Bild 7). FD-Verfahren mit TVD-Eigenschaften sind inzwischen Stand der Technik. Allerdings hat die Sache einen Haken: Wie Jonathan Goodman vom Courant-Institut und Randall LeVeque von der Universität Seattle (US-Bundesstaat Washington) 1985 zeigten, muß man in mehr als einer Raumdimension für die TVD-Eigenschaft so viel Genauigkeitsverlust in Kauf nehmen, daß man auf den Stand der einfachsten Differenzenverfahren aus von Neumanns Zeiten zurückgeworfen wird.

Es ist ein (schwacher) Trost, daß auf diesem Feld die Numerik nicht allzuweit hinter der Theorie zurückbleibt: Es gibt ohnehin kaum beweisbare Aussagen über das Verhalten der Lösung mehrdimensionaler Probleme. Darum treiben wir das operator splitting noch einen Schritt weiter: Wir behandeln die Anteile des Transports in den drei Raumrichtungen nacheinander und unabhängig voneinander – also jeweils eindimensional –, so wie schon Transport und Diffusion getrennt behandelt werden.

Da man den Widerspruch zwischen den Forderungen nach Genauigkeit und TVD-Eigenschaft bereits Mitte der achziger Jahre erkannt hatte, tauchten bald darauf neue Verfahren auf, die der Totalvariation der numerischen Lösung erlaubten, in gewissen kleinen Grenzen im Laufe der Zeit zu wachsen. Diese ENO-Verfahren (essentially non-oscillatory) ermöglichen die Konstruktion von – im Prinzip – beliebig genauen Methoden, sind aber in der Anwendung aufwendig und deshalb in der Praxis noch zu teuer.


Lösung von Diffusions- und Reaktionsgleichung

Im Vergleich zu dem großen theoretischen Apparat, der für die Transportgleichungen existiert, gibt es zu den Diffusions- und den Reaktionsanteilen relativ wenig zu sagen. Diffusionsterme sind numerisch eher gutmütig, erfordern zu ihrer Diskretisierung keinen besonderen Scharfsinn und glätten die Lösung – es sei denn, Turbulenz kommt ins Spiel.

Für dieses Phänomen ist die theoretische Basis noch äußerst dünn. Eigentlich kann man einem Turbulenzmodell nur dann trauen, wenn man Rechenergebnisse mit Experimenten vergleichen kann. Andererseits enthalten detaillierte Modelle zusätzliche partielle Differentialgleichungen, so daß sich der Aufwand stark erhöht; häufig treten numerische Instabilitäten auf. Die dadurch verursachte Unsicherheit beeinträchtigt die Qualität der errechneten Lösung im ganzen. Bessere Turbulenzmodelle sind dringend erforderlich; Fortschritte erwarten wir von direkten Simulationen des Turbulenzvorgangs.

Bei den chemischen Reaktionen liegt, wie erwähnt, die wesentliche Schwierigkeit in der Steifigkeit der zugehörigen Differentialgleichungen, die durch sehr unterschiedliche Reaktionskonstanten verursacht wird. Entgegen früheren Befürchtungen ist dieses Problem jedoch längst nicht so ausgeprägt wie etwa bei Verbrennungsprozessen im Kolbenmotor; klassische Verfahren zur Lösung reichen in den meisten Fällen aus. Zudem gibt es für steife Systeme aus der jüngeren Vergangenheit eine gut ausgebaute Theorie mit zugehörigen Verfahren.

Chemische Reaktionen sind lokal: Was sich an einem Punkt des Rechengebietes abspielt, hängt nur vom Zustand in diesem Punkt ab und nicht von dessen Umgebung. Aus diesem Grunde ist ein Verfahren anwendbar, das für nicht-lokale Terme in mehreren Raumdimensionen einen extrem hohen Rechenaufwand erfordern würde: die implizite Diskretisierung bezüglich der Zeit.

Diskretisieren ist das Annähern einer Differentialgleichung durch eine Differenzengleichung; deren Lösung wiederum ergibt den (genäherten) Zustand des Systems für den neuen Zeitpunkt. Stellt man die Differenzengleichung für den alten Zeitpunkt auf, ist das Lösen sehr einfach. Wenn man sie dagegen für den neuen Zeitpunkt aufstellt, geht der neue, noch unbekannte Zustand an mehreren Stellen in die Gleichung ein; die Information ist nur noch implizit in der Gleichung enthalten und muß mit entsprechendem Aufwand aus ihr extrahiert werden (Bild 4 unten). Deshalb heißt das zugehörige Verfahren implizit im Gegensatz zu den einfacheren expliziten Verfahren. Andererseits hat dieses mehrfache Einbringen der Unbekannten einen stabilisierenden Effekt: Implizite Verfahren erlauben wegen ihrer höheren Stabilität größere Zeitschritte als explizite.

Im Prinzip muß man die zeitliche Schrittweite so wählen, daß der am schnellsten ablaufende Prozeß noch mit ausreichender Genauigkeit erfaßt wird. Indem man die Reaktionsterme implizit in der Zeit diskretisiert, erreicht man, daß der Zeitschritt nicht allzu klein und der Rechenaufwand pro Zeitschritt nicht allzu groß wird.


Anwendungen

Bei der DLR wurden zahlreiche der dargestellten einzelnen Ideen in zwei Programmpaketen realisiert, die auf regelmäßigen beziehungsweise unregelmäßigen Gittern basieren. Der t-Code (der griechiche Buchstabe tau steht für triangular, adaptive, upwind) ist ein FV-Verfahren auf unregelmäßigen Dreiecksnetzen; es fügt je nach Bedarf Punkte ins Netz ein oder eliminiert sie. Es verwendet einen approximativen Riemann-Löser. Zusätzliche Genauigkeit gewinnt er durch einen Nachbearbeitungsschritt, der im Endeffekt darauf hinausläuft, daß nicht mit stückweis konstanten, sondern mit stückweis linearen Funktionen gearbeitet wird – die gleichwohl an den Grenzen der (drei- oder viereckigen) Rechenzellen unstetig sein dürfen.

Es gibt Standardprobleme, an denen man die Qualität eines neuentwickelten Codes zu messen pflegt. Für Überschallströmungen ist eines von ihnen die reibungsfreie Strömung in einem Kanal mit einspringender Stufe (forward-facing step). Die Aufgabe ist so ausgelegt, daß eine der drei Raumrichtungen, nämlich horizontal quer zur Strömungsrichtung, keine Rolle spielt; die Reduktion auf zwei Raumdimensionen verkürzt die Rechenzeit erheblich. Bild 8 zeigt die Lösung durch den t-Code.

Inzwischen können wir routinemäßig reagierende Über- und Hyperschallströmungen untersuchen, die für moderne Raumflugsysteme eine Rolle spielen.

Das Übungsbeispiel zeigt bereits wesentliche Eigenschaften von Hyperschallströmungen: Es bildet sich zunächst eine Bugstoßwelle; weitere Stoßfronten entstehen dort, wo der Rand des umströmten Körpers nicht glatt ist, sondern einen Knick hat – bei einem Überschallflugkörper dort, wo Flügel oder Steuerklappen am Rumpf ansetzen oder wo die Pilotenkanzel mitsamt den Sichtfenstern aus dem ansonsten stromlinienförmigen Rumpf hervorspringt.

Stoßfronten werden in komplizierter Weise an festen Wänden und zum Teil auch aneinander reflektiert. Wo eine Stoßfront auf die Oberfläche des Flugkörpers trifft, entstehen starke Hitze und großer Druck mit zerstörerischer Wirkung. Die Form ist also so zu gestalten, daß diese Zustände in keiner Fluglage auftreten. Zudem möchte man nach den Erfahrungen beim ersten Raumfähren-Flug verständlicherweise sicherstellen, daß die Steuerklappen wirksam sind.

Um erste Erfahrungen mit (noch zu entwerfenden) numerischen Verfahren zu gewinnen, wurde im Rahmen des Programms Hermes eine vereinfachte Modellkonfiguration als Übungsaufgabe definiert. Sie besteht aus einem Hyperboloid mit angesetzter Rampe (hyperboloid flare; Bild 9). Damit kann man die Wechselwirkungen zwischen Verdichtungsstößen und der Grenzschichtströmung im Bereich der Steuerklappe unter Wiedereintrittsbedingungen studieren, das heißt bei hohen Machzahlen, großen Anstellwinkeln und unter Berücksichtigung von Hochtemperatureffekten.

Wie um den vollständigen Flugkörper bildet sich auch um den gesamten Vorkörper des hyperboloid flare ein Verdichtungsstoß, der mit dem lokalen Stoßsystem, das durch die Umlenkung der Strömung an der Rampe hervorgerufen wird, in Wechselwirkung tritt. Dadurch entsteht ein Gebiet mit abgelöster Strömung, das heißt, die Stromlinien verlaufen nicht parallel zum Körper. Es zeigt sich, daß die Strömungsgrößen unter verschiedenen Annahmen für die chemischen Reaktionen völlig verschiedene Werte annehmen (Bild 9 unten). Insbesondere führen die früher üblichen vereinfachenden Annahmen einer gefrorenen beziehungsweise einer im chemischen Gleichgewicht befindlichen Strömung grob in die Irre – auch, was die Wirksamkeit einer Steuerklappe angeht.

Der Bugstoß kann stromab in Wechselwirkung mit dem abgelösten Stoß vor einem Flügel oder Leitflügel treten. Wieder will man vermeiden, daß davon eine zerstörerische Wirkung auf die Oberfläche ausgeht. Das vereinfachte Testproblem ist für diesen Fall ein quer angeströmter Zylinder, auf dessen Bugstoß ein weiterer schräger Verdichtungsstoß trifft (Bild 10).

Inzwischen ist ein bei der DLR in Göttingen entwickeltes Pilotverfahren für reagierende Strömungen mit einem bei der DLR in Braunschweig entwickelten FV-Verfahren kombiniert worden, das durch Verwendung von Mehrgitter-Prinzipien äußerst effizient ist (vergleiche Spektrum der Wissenschaft, April 1990, Seite 78). Das Programm CEVCATS-N (CEll Vertex Central Averaging Time-stepping Scheme – Nonequilibrium) berechnet nicht nur zweidimensionale Modellströmungen, sondern die volle dreidimensionale Strömung um ein Wiedereintritts-Fluggerät – einschließlich des eingangs erwähnten Problems mit der Steuerklappe (Bild 1).

Die beschriebenen Algorithmen sind nicht nur für die ursprünglich anvisierten Zwecke geeignet, sondern zum Beispiel auch für Verbrennungsprozesse: vom Rasenmähermotor bis zum Triebwerk eines in der Entwicklung befindlichen Überschall-Verkehrsflugzeugs. Weitere Anwendungsfelder sind moderne Wasserstoff-Antriebssysteme für schubstärkere Versionen der Trägerrakete Ariane, die Satelliten zu mehreren gebündelt und damit kostengünstiger in ihre Umlaufbahn bringen sollen, neue Verkehrsflugzeuge mit geringerem Treibstoffverbrauch und die aerodynamische Optimierung künftiger Hochgeschwindigkeitszüge. Die Navier-Stokes-Gleichungen erfassen auch diese Phänomene. Prinzipiell ist den numerischen Verfahren nichts hinzuzufügen – im Detail dagegen allerlei.

Bislang werden bei der Behandlung der reagierenden Strömungen nur vereinfachte Reaktionsschemata verwendet, weil die Gesamtheit der im Prinzip möglichen chemischen Reaktionen noch jeden Hochleistungscomputer überfordert. Trotzdem muß das globale Verhalten der ablaufenden Reaktionen richtig beschrieben werden. Hier ist noch viel zu tun.

Für die Turbulenz fehlt bislang ein breit anwendbares physikalisches Modell, das in praktikable numerische Verfahren umsetzbar wäre. Da indes die meisten technisch interessierenden Verbrennungsvorgänge auf einer turbulenten Flamme beruhen, gilt es auch dafür geeignete Modelle zu entwickeln, in das numerische Verfahren zu implementieren und anhand von Experimenten zu validieren. Gleiches gilt für die Modellierung der Tröpfchenbildung, die berücksichtigt werden muß, wenn sich flüssige und gasförmige Treibstoffe vermischen. Auf diesem Gebiet der Modellbildung werden derzeit in allen erwähnten Bereichen große Anstrengungen unternommen. Der Erfolg der numerischen Berechnung von komplexen Strömungen wird nicht zuletzt von der Qualität der verwendeten Modelle abhängen.

Literaturhinweise

- Hypersonic and High Temperature Gas Dynamics. Von John D. Anderson jr. McGraw-Hill, 1989.

– Non-Equilibrium Hypersonic Aerothermodyamcs. Von Chul Park. Wiley, 1990.

– Numerical Simulation of Shock/Shock and Shock-Wave/Boundary-Layer Interactions in Hypersonic Flows. Von G. Brenner, T. Gerhold, K. Hannemann und D. Rues in: Computers and Fluids, Band 22, Heft 4/5, Seiten 427 bis 439, 1993.

– On the Design of an Upwind Scheme for Compressible Flow on General Triangulations. Von Th. Sonar in: Numerical Algorithms, Band 4, Seiten 135 bis 149, 1993.

– Experimental Verification of Real Gas Effects in High-Enthalpy Flows. Von Georg Eitelberg in: Shock Waves @ Marseille. Proceedings of the 19th International Symposium on Shock Waves, Marseille, Juli 1993. Springer, Berlin 1995.

– Efficient Numerical Simulation of Complex 3D Flows with Large Contrast. Von R. Radespiel, J. M. A. Longo, S. Brueck und D. Schwamborn in: 77th AGARD Fluid Dynamics Panel Meeting and Symposium on Progress and Changes in CFD Methods and Algorithms, Sevilla, Oktober 1995.


Aus: Spektrum der Wissenschaft 7 / 1996, Seite 72
© Spektrum der Wissenschaft Verlagsgesellschaft mbH