Numerische Simulation – von der Formel zum bunten Bild

Am Schweizerischen Institut für Informationswissenschaft (SII) der Fachhochschule Graubünden (FHGR) wurde mit DAViS ein Kompetenzzentrum gegründet, das die Themenschwerpunkte Datenanalyse, Visualisierung und Simulation vereint und zusammen mit Partnern aus Forschung und Industrie sowie dem Swiss National Supercomputing Centre (CSCS) in Lugano gemeinschaftliche Projekte zu aktuellen Forschungsfragen bearbeitet. Exemplarisch soll hier der Themenschwerpunkt Simulation herausgegriffen und ein Projekt aus dem Bereich der Hochwasserereignisse vorgestellt werden.

Mit der gestiegenen Rechenleistung moderner Computer im Lauf der letzten Jahrzehnte hat sich die numerische Simulation neben Theorie und Experiment als „dritte Säule“ des Erkenntniserwerbs fest etabliert. Dabei werden mithilfe von mathematischen Modellen unterschiedliche Szenarien am Rechner getestet (simuliert), um quantitative Vorhersagen treffen zu können. Typische Vertreter numerischer Simulation sind Strömungssimulationen etwa zur Bestimmung des Wasserpegels bei Überschwemmungen, für Klimavorhersagen oder die Umströmung einer Tragfläche zur Optimierung von Auftrieb und Luftwiderstand sowie Festkörpersimulationen etwa zur Strukturanalyse von Bauwerken, für Crashtests oder im Rahmen biomechanischer Anwendungen zur Vorhersage von Materialverschleiß an Knochen und Wirbelkörpern. Die numerische Simulation ist immer dann vorteilhaft, wenn Experimente zu aufwendig sind (d. h. zu teuer und/oder zu zeitintensiv), nicht durchgeführt werden können (z. B. Supernovae) oder nicht durchgeführt werden sollen (z. B. Lawinenabgänge).

Grundlage der numerischen Simulation sind in der Regel mathematische Modelle, die physikalische Zusammenhänge mithilfe von Differentialgleichungen beschreiben. In der Strömungsmechanik sind bspw. die Euler-Gleichungen sowie deren Erweiterung auf linear-viskose newtonsche Flüssigkeiten und Gase weit verbreitet; letztere sind auf die beiden Mathematiker und Physiker Claude Louis Marie Henri Navier und George Gabriel Stokes zurückzuführen und werden daher als Navier-Stokes-Gleichungen oder kurz NS-Gleichungen bezeichnet. Sie bestehen aus einer Kontinuitätsgleichung (i) zur Gewährleistung der Divergenzfreiheit (d. h. Quellfreiheit des Rechengebiets) sowie aus einer Impulsgleichung (ii) je Raumdimension. Vereinfacht können die NS-Gleichungen für inkompressible Fluide nach Hirsch 2007 folgendermaßen aufgeschrieben werden,

⋅u→=0,(i)

∂∂t+(u→⋅)u→=νu→-1ρp+g→,(ii)

wobei u→ die Geschwindigkeit, p den Druck, ρ die Dichte, ν die kinematische Viskosität und g→ die Erdbeschleunigung ausdrücken. Je nach gewählter Viskosität lassen sich damit zähflüssige Fluide (z. B. Öl oder Honig) bis hin zu Luft und Wasser abbilden.

Da für die NS-Gleichungen keine analytischen Lösungen existieren, müssen diese numerisch gelöst werden. Hierzu ist es jedoch erforderlich, die Gleichungen per se sowie das Rechengebiet, innerhalb dessen mithilfe der Gleichungen eine Lösung des zu Grunde liegenden Strömungsproblems gefunden werden soll, zu diskretisieren. Anstelle der exakten Lösung in unendlich vielen Punkten werden die Gleichungen nur an endlich vielen Punkten ausgewertet und damit eine diskrete Näherungslösung bestimmt. Infolgedessen ergibt sich ein (Diskretisierungs)Fehler, der die Qualität der Lösung maßgeblich beeinflusst und von der Dichte der Punkte, genauer vom Abstand zwischen zwei Punkten abhängt. Allgemein lässt sich festhalten, je kleiner der Abstand, desto besser die Qualität, aber auch desto höher der Rechenaufwand. Zur Diskretisierung der Gleichungen stellt die Mathematik ein reichhaltiges Angebot an Werkzeugen zur Verfügung, auf die an dieser Stelle nicht weiter eingegangen werden soll. Der interessierte Leser wird verwiesen an Ferziger und Perić (2008) oder Knabner und Angermann (2000).

Analog zu den Gleichungen muss das Rechengebiet mit der zu umströmenden Geometrie diskretisiert, d. h. in endlich viele Teilstücke zerlegt werden. Auch hier bietet der Methodenbaukasten eine Vielzahl an unterschiedlichen Verfahren (Samet 1990), stellvertretend soll hier das raumpartitionierende Verfahren der Spacetrees (Mundani 2006) vorgestellt werden. Ausgehend von einem Rechteck im 2D oder einem Quader im 3D wird dieses Gebiet – zusammen mit der darin enthaltenen Geometrie – durch rekursive Bisektion entlang jeder Raumrichtung so lange in kongruente Teilstücke zerlegt, bis die resultierenden Elemente (im Englischen spricht man von Volume Element, abgekürzt Voxel) komplett innerhalb oder außerhalb der Geometrie liegen. Abbildung 1 zeigt die ersten Verfeinerungsschritte einer einfachen 2D-Geometrie, Abbildung 2 das Ergebnis der Diskretisierung einer Kraftwerksgeometrie bei unterschiedlichen Auflösungen.

Abbildung 1

Ausgangssituation (links) im 2D gefolgt von den ersten drei Verfeinerungsschritten.

Abbildung 2

Kraftwerksgeometrie^[1] (links) nach acht (Mitte) und zehn (rechts) Verfeinerungsschritten.

Dabei haben Spacetrees deutliche Vorteile gegenüber klassischen Verfahren wie dem Normzellenaufzählungsschema, bei dem das Rechengebiet in äquidistante Zellen – also mit gleicher Seitenlänge – zerlegt wird. Für eine Seitenlänge h ergibt das n=1/h Zellen je Raumrichtung, was folglich für das Normzellenaufzählungsschema im 2D und 3D zu einer Komplexität von O(n²) bzw. O(n³) Zellen führt. Da Spacetrees nur entlang der Oberfläche einer Geometrie verfeinern, führt dies bei gleicher geometrischer Genauigkeit (d. h. die feinsten Zellen haben die gleiche Größe) zu einer verbesserten Komplexität von O(n) und O(n²) im 2D bzw. 3D – fraktale Geometrien einmal außer Acht gelassen.^[1]

Nach erfolgter Diskretisierung des Rechengebiets muss für jede Zelle noch festgelegt werden, welche Eigenschaften diese besitzt. Typischerweise wird hier zwischen Fluid- und Hinderniszellen unterschieden; auf ersteren werden die NS-Gleichungen ausgewertet, um eine Näherungslösung zu bestimmen, während auf letzteren die Lösung in Form von Randbedingungen vorgegeben ist. Zur Bestimmung der Näherungslösung bedarf es nunmehr eines numerischen Verfahrens, das – im Fall transienter, d. h. instationärer Probleme – zu jedem (diskreten) Zeitschritt für jede Zelle die gesuchten Lösungen für Druck p und Geschwindigkeit u→ berechnet.

Durch Umstellung von Gleichung (ii) mithilfe der Chorin’schen Projektionsmethode (Chorin 1967) erhält man schließlich zwei Gleichungssysteme für Druck und Geschwindigkeit, die sich nacheinander lösen lassen. Dabei werden zur Berechnung des Drucks p^t+1 im Zeitschritt t+1 nur die bereits bekannten ‚alten‘ Geschwindigkeiten u→t vom vorherigen Zeitschritt t benötigt, bevor mittels ‚neuem‘ Druck anschließend die Geschwindigkeiten aktualisiert werden und man u→t+1 erhält. Die Lösung der Druckgleichung ist dabei um Größenordnungen aufwendiger und bedingt neben effizienten numerischen Verfahren zumeist auch Konzepte des Parallelrechnens, bei dem das zu Grunde liegende Gleichungssystem auf mehrere Rechner verteilt gelöst wird. Auch hier bietet die Mathematik einen ganzen Zoo von Methoden an, die von direkten Lösern wie der LR-, Cholesky- oder QR-Zerlegung über iterative Verfahren wie Jacobi, Gauß-Seidel, CG oder Mehrgitter bis hin zu Spektralmethoden reichen. Nicht jede Methode ist für jedes Problem geeignet, manche lassen sich einfacher umsetzen und manche besser parallelisieren. Weiterführende Literatur findet der interessierte Leser etwa unter Meister (2015) sowie Golub und Van Loan (2013).

Wie bereits oben aufgeführt, lassen sich manche Problemstellungen nur durch den Einsatz von Großrechnern oder Supercomputern lösen. Ausschlaggebend für die Rechenkomplexität ist die gewählte Auflösung des Rechengebiets, d. h. die Anzahl der Zellen. Bereits bei einer äquidistanten Diskretisierung von 100×100×100 Zellen ergibt das für den Druck ein Gleichungssystem mit n=10⁶ Unbekannten. Löst man dieses mit der sicherlich nicht kompetitiven Schulbuchmethode (LR-Zerlegung oder umgangssprachlich Gauß’sche Elimination), so führt das zu einer Rechenkomplexität von O(n³) Operationen, für die ein Standardcomputer ohne Parallelisierung im – unrealistischen – Idealfall etwa 10 Jahre benötigt.

Aktuell (Stand November 2019^[2]) erreicht der schnellste Rechner der Welt (Summit am Oak Ridge National Laboratory, USA) auf insgesamt 2,4 Millionen Kernen eine Rechenleistung von 148,6 PFlop/s – das sind 148.600.000.000.000.000 Gleitkommaoperationen pro Sekunde. Der schnellste Rechner der Schweiz, der Piz Daint am Swiss National Supercomputing Centre, kommt immerhin auf 21,2 PFlop/s bei knapp 400.000 Kernen. Die Lösung des vorherigen Problems mit n=10⁶ Zellen würde auf dem Summit bei gleichem Lösungsverfahren im – unrealistischen – Idealfall etwa 7 Sekunden benötigen. Mit unserem hauseigenen Strömungscode MPFluid konnten wir auf einem der drei Höchstleistungsrechner Deutschlands (SuperMUC am Leibniz-Rechenzentrum in Garching bei München) auf gut 140.000 Kernen ein Strömungsproblem mit 700 Milliarden Unbekannten rechnen.

Dabei stand eine Hochwassersimulation im Fokus, wie sie etwa nach Starkregenereignissen auftreten kann. Im vorliegenden Fall wurde die Münchner Innenstadt virtuell geflutet, um im Anschluss aufgrund vorhergesagter Pegelstände eine konkrete Schadensbemessung durchführen zu können. Grundlage für die Simulation war ein hoch detailliertes BIM-Modell (Building Information Modelling) (Eastman, Teicholz, Sacks und Liston 2011) der Technischen Universität München, das in ein GIS-Modell (Geographische Informationssysteme) (Brinkhoff 2013) der Innenstadt rund um den Hauptbahnhof eingebettet wurde. Über das GIS-Modell waren zudem Angaben über unterirdische Bauwerke wie etwa die Kanalisation oder Regenrückhaltebecken vorhanden. Unter der Annahme, dass nach einem Starkregen die Rückhaltebecken sowie die Kanalisation komplett mit Wasser gefüllt waren, wurde eine mehrskalige Simulation auf unterschiedlichen Auflösungsebenen durchgeführt, d. h. die Ergebnisse einer gröberen Ebene fungierten als Randbedingung für eine feinere Ebene. Auf der gröbsten Ebene (Stadtlevel) wurden dabei globale Effekte, d. h. die Ausbreitung der Flutwelle durch die Straßen betrachtet, auf der mittleren Ebene (Gebäudelevel) konnten lokale Effekte, d. h. das Eindringen der Flutwelle in einzelne Gebäude berücksichtigt werden, und auf der feinsten Ebene (Detaillevel) ließen sich schließlich konkrete Effekte wie die Höhe des Wasserpegels im Innenraum quantifizieren.

Die folgenden Abbildungen 3 bis 6 zeigen Ergebnisse der Hochwassersimulation im Münchner Innenraum zu unterschiedlichen Zeitschritten auf unterschiedlichen Auflösungen. Die markierten Bereiche stellen die jeweils feiner aufgelösten Regionen auf der nächsten Detailebene dar.

Abbildung 3

GIS-Modell der Münchner Innenstadt mit darunter liegender Kanalisation.

Abbildung 4

Bestimmung globaler Effekte (Stadtlevel) am GIS-Modell.

Abbildung 5

Berücksichtigung lokaler Effekte (Gebäudelevel) am BIM-Modell.

Abbildung 6

Quantifizierung konkreter Effekte (Detaillevel) zur Schadensvorhersage.