Werkzeug

Da diese Web-Seite für Modellflieger entstanden ist, werfe ich jetzt nicht mit Partiellen Differentialgleichungen um mich, sondern versuche verbal zu erklären, welche Variablen berechnet werden. Die Gleichungen teilen sich auf in einen Teil für die mittlere Strömung und in einen für die Turbulenzschliessung. Eigentlich ist es ganz einfach: Man formuliert die Erhaltungssätze für Impuls, Masse und Energie inkl. Newton's Kraft gleich Masse mal Beschleunigung und eventuell garniert mit einer Viskosität für Luft und errechnet eine Verteilung von Geschwindigkeiten und Druck. Die sog.  Turblenzschließung ist dann verantwortlich, via Vermischung eine Grenzschicht zu erzeugen und verantwortlich ist auch für den Luftwiderstand und den Auftrieb. 

Bewegungsgleichung

Die Bewegungsgleichungen für die mittlere Strömung werden durch 3 prognostische Gleichungen für den Impuls (Dichte mal Geschwindigkeit), der Dichte über eine sog. Kontinuitätsgleichung, der potentiellen Temperatur, einer diagnostischen Gleichung genannt Zustandsgleichung, die die Beziehung zwischen Druck, Temperatur und Dichte herstellt und einer Gleichung, mit der man potentielle Temperatur in "in situ" Temperatur umrechnen kann, numerisch gelöst. Komplettiert werden die 3 Gleichungen für den Impuls neben den Termen für Impulstransport, Druckgradienten und Schwerkraft mit einer Impulsdiffusion, die je nach Richtung mit räumlich variablen Diffusionskoeffizienten in der Lage sind, die Auswirkungen der Turbulenz abzubilden. Turbulenz ist die Voraussetzung für die Bildung einer Grenzschicht, d.h. eines Geschwindigkeitsprofils an einer festen Oberfläche, und davon bietet ein Flugzeug genug. Die Gleichungen für die mittlere Strömung zusammen mit den dazugehörigen Parametern ergeben ein in sich geschlossenen Gleichungssatz.

Turbulenzschließung

Die Turbulenzschliessung selbst bestand ursprünglich aus 3 Gleichungen, nämlich für horizontale und vertikale turbulente kinetische Energie und für die Längenskala der erzeugten Wirbel. Ich bezeichne dies als "2k-L"-Modell. Auf internationaler Ebene findet man diverse Ansätze, z.B. "k-kL" neben "k-w", "k-ε," und neben anderen meist einfacheren. Grundsätzlich muss man feststellen, dass es sich bei den diversen Varianten von Turbulenzschliessungen nicht um eine geschlossene Theorie handelt. Stattdessen zeigen die verbleibenden Konstanten, die aus Experimenten zu bestimmen sind, dass es sich um Näherungen handelt, wenn auch mit hohem theoretischem Geschick erreicht. Meine Variante der Turbulenzschliessung geht zurück auf einen ehemaligen Kollegen Lautenschlager (1983 Techn. Report 83/E/60, GKSS, Geesthacht; 1985, Techn. Report 85/E/18, GKSS, Geesthacht), der diese Form der Turbulenzschliessung für großräumige Anwendungen eingeführt hat. Die 3 Gleichungen beinhalteten den Transport und die Entstehung von Turbulenz aufgrund der Deformation des Strömungsfeldes, der Dissipation und Umverteilung zwischen der horizontalen und vertikalen Komponente der turbulenten, kinetischen Energie und einen Term anschaulich ausgedrückt für die Alterung der Turbulenz, die als Änderung der Längenskala der Turbulenz mit der Zeit abgebildet wird. Erstaunlich für mich war, dass die 4 experimentell zu bestimmenden Konstanten bis auf eine bereits 1985 nahezu die Werte hatten, die heute nach Umrechnung auf das "k-kL"-Konzept von Khaled S. Abdol-Hamed (NASA Langley Research Center, Hampton VA, USA) in einem Artikel in "International Journal of Aerospace Engeneering, Article ID 987682") 2015 veröffentlicht wurden. Die vierte Konstante war vergleichsweise robust zu bestimmen und ist zu verstehen als Anpassung an kleinskalige Strömungen - eine Skala von Millimetern anstatt Kilometern bedeutet, dass die Brunt-Väisäla-frequency hier nichts zu suchen hat - weil vernachlässigbar. Nachdem mir bei einigen Experimenten die Ergebnisse der Turbulenzschließung missfallen haben, habe ich in Teilen das Vorgehen von Khaled S. Abdol-Hamed (NASA Langley Research Center, Hampton VA, USA) berücksichtigt. Aber, nachdem ich auch mit diesen Ergebnissen nicht einverstanden war, habe ich ein "3K-L"-Konzept erstellt - mit vielen weiteren - neutral formuliert - Änderungen. Die Gleichungen für die Turbulenzschliessung ist zu finden hier zusammen mit den dazugehörigen Parametern. Um Fragen zu relativieren, warum ich kein K-ε-Konzept verwende, verweise ich auf folgende Analyse, die belegt, dass K-ε, K-KL und K-L mathemtisch identisch sind, solange man die jeweiligen Konstanten aus dem einen Konzept  in ein anderes umrechnet.

Geschichte des CFD-Modells

Der Kern des Modells geht zurück auf Oberhuber et al. (Volcanic plume simulation on large scales, Journal of volcanology and geothermal research, 1998, 29-53) und wurde von Doktoranden ergänzt durch Parametrisierungen der Mikrophysik, genauer Wolkenphysik, unter Beteiligung von Asche und allen Aggregatszuständen von Wasser mit dem Ziel, realistisch die Entwicklung einer Aschewolke bei einem Vulkanausbruch zu simulieren. Voraussetzung dazu ist, dass alle transportierten Stoffe eine beliebig hohe Konzentration besitzen dürfen und somit dynamisch aktiv sind. Dieses Modell ist meine Wissensbasis für die auf dieser Internetseite verwendeten Ergebnisse meines privat entwickelten und ausschließlich zu privaten Zwecken verwendeten virtuellen Windkanals. Eigentum obigen Modells ist das Max-Planck-Institut für Meteorologie und die Deutsche Klima-Rechenzentrum GmbH, beide in Hamburg, und hat mit meiner Version des virtuellen Windkanals inzwischen wenig gemeinsam. So wurden die offenen Randbedingungen neu konzipiert, alle Diskretisierungen der Transportterme so ausgetauscht, dass Impuls, Masse und Wärme exakt erhalten werden mit dem Ziel auch quantitativ der Realität näher zu kommen oder z.B. die Aufteilung der turbulenten kinetischen Energie auf die Verhältnisse an einem Modellflugzeug angepasst, und grundlegende Änderungen bzgl. Programmstruktur, Performancemaßnahmen, Konvergenz der iterativ zu lösenden implizit in der Zeit formulierten Gleichungen verbessert und Anpassungen durchgeführt, die ein multi-threading basierend auf openmp ermöglichen. 

Numerische Umsetzung

Die Lösungsmethode für ein implizites Zeitschrittverfahren basiert auf der prinzipiellen Idee von Kwizak und Robert (1971), aus den Gleichungen für Geschwindigkeit und Druck eine Wellengleichung abzuleiten, die dann gelöst wird, um schließlich aus dem sich ergebenden Druck zum neuen Zeitpunkt die neuen Geschwindigkeiten zu ermitteln. Der Pferdefuß ist allerdings, dass dies nur für einen Gleichungstyp wie Flachwassergleichung einfach unmsetzbar ist. Hier muss etwas mehr struktureller Aufwand getrieben werden, um letztendlich für alle 9 prognostischen Variablen ein implizites Zeitschrittverfahren umzusetzen. Knackpunkt ist dann aber die iterative Lösung für x-Millionen Gitterpunkte, die massiv durch den Einsatz von linearen Gleichungssystemen unterstützt wird. Ziel ist es dann, mit hohen CFL-Zahlen eine stabile Lösung zu ermitteln. Für Geschwindigkeiten verwende ich situationsgesteuert eine maximale CFL-Zahl von bis zu 3. Schallwellen breiten sich dann typisch mit CFL=30 aus. Die Gleichungen werden inkl. der Randbedingungen auf einem Arakawa C-Gitter diskretisiert, veranschaulicht durch folgende Darstellung. Vorteil ist, dass die algebraischen Gleichungen, die eine Näherung an die nicht lösbaren, partiellen Differentialgleichungen sind, keine Filter benötigen, um die Lösung stabil zu halten.