Autorenliste · 2012. 6. 6. · 4.3.1 Validierung der Monte Carlo Simulation der Partikeldynamik...

Autorenliste:

Forschungsstelle 1:

E. Kruis

B. Sherif

J. Wei

Numrax:

R. Vilsmeier

Forschungsstelle 2:

S. Haep

T. Zeiner

T. van der Zwaag

Forschungsstelle 3:

H. Orthner

C. Weise

H. Wiggers

I. Wlokas

Schlussbericht 09/2009 – 01/2012 Kombination CFD – Monte Carlo i

Inhalt:

1 Wissenschaftliche Problemstellung 1

1.1 Ausgangssituation 1

1.2 Anlass für den Forschungsantrag 2

2 Forschungsziel 3

3 Zusammenfassung der Ergebnisse 3

4 Lösungsweg 6

4.1 Modellentwicklung 6

4.1.1 Grundlegende Problemstellung 6

4.1.1.1 Präambel zum Rechenleistungsbedarf und numerischen Aufwand 7

4.1.1.2 Begründung für die Benützung von Grafikkarten für die Simulation von praxisnahen Reaktoren 9

4.1.1.3 Vorgehensweise 10

4.1.2 Optimierung der Monte Carlo Simulationstechnik für die Partikeldynamik 12

4.1.3 Partikeltransportmodul für die Monte Carlo-Simulation auf Gittern beliebiger Struktur 21

4.1.3.1 Eulerscher versus Lagrangescher Partikeltransport 21

4.1.3.2 Eulerscher Partikeltransport 22

4.1.3.3 Lagrangescher Partikeltransport 23

4.1.3.4 Bestimmung der Zielzellen 26

4.1.3.5 Algorithmische Struktur und Auslegung der Rechenanlagen 28

4.1.4 Anpassung der Monte Carlo Software zur Beschreibung willkürlicher Partikelprodukteigenschaften 37

4.1.5 Konzeption und Realisierung eines preiswerten Hochleistungsrechnersystems optimiert für Monte Carlo-basierte Partikeldynamik und -transport 39

4.1.6 Deterministische Modelle 48

4.1.6.1 Deterministische Beschreibung der Fluiddynamik 48

4.1.6.2 OpenFOAM 51

4.1.6.3 Programmerweiterung zur Beschreibung der Partikelbildung 51

Schlussbericht 09/2009 – 01/2012 Kombination CFD – Monte Carlo ii

4.2 Experimente zur Nanopartikelsynthese 55

4.2.1 Wissenschaftliche Problemstellung 55

4.2.2 Forschungsziel 55

4.2.3 Lösungsweg 56

4.2.3.1 Aufbau des Heißwand-Reaktors 56

4.2.3.2 Auswahl der Betriebsparameter 58

4.2.3.3 Partikelprobenahme und Messung 60

4.2.3.4 Pyrolyse von Fe(CO)5 und Bildung von Clustern 64

4.3 Validierung der Methoden und Modelle 68

4.3.1 Validierung der Monte Carlo Simulation der Partikeldynamik auf der GPU aufgrund Koagulation 68

4.3.2 Validierung der Monte Carlo Simulation von gleichzeitigem Transport und Partikeldynamik 70

4.3.3 Validierung der CFD-Methode am Beispiel eines virtuellen Testreaktors 74

4.4 Anwendung auf praxisnahe Fallbeispiele 78

4.4.1 CFD-Simulationen für den Heißwandreaktor von FST 3 78

4.4.2 Validierung des CFD-Simulationsverfahrens 79

4.4.3 Variation der Prekursorkonzentrationen 88

4.4.4 Kombinierte CFD- Monte Carlo Simulation des Heißwandreaktors aus FST 3 90

5 Gegenüberstellung der Ergebnisse mit den Zielsetzungen des ursprünglichen Forschungsantrages 92

6 Darstellung des wissenschaftlich-technischen und wirtschaftlichen Nutzens der erzielten Ergebnisse insebesondere für KMU 97

7 Ergebnistransfer in die Wirtschaft und Liste der Veröffentlichungen 98

8 Durchführende Forschungsstellen 101

9 Literatur 102

Schlussbericht 09/2009 – 01/2012 Kombination CFD – Monte Carlo 1

1 Wissenschaftliche Problemstellung

1.1 Ausgangssituation

Die großtechnische Herstellung von Nano- und Mikro-Partikeln, wie z. B. die Pigment- und Rußproduktion, erfolgt überwiegend in Aerosolreaktoren und Kristallisatoren. Die Eigenschaften der dispersen Phase im Reaktor sind aus zwei Gründen schwer vor-auszusagen: 1) Die Zerfalls-, Bildungs- und Wachstums- bzw. Agglomerationsprozes-se in den Synthesereaktionen überlagern sich, so dass die Mischung und Transport der Reaktanten sowie die Kinetik der Prozesse nur über ort- und zeitaufgelöste de-terministische Modelle simuliert werden kann. 2) Bisher ist kommerzielle Software (bspw. Fluent®) nur in der Lage, eine Partikeleigenschaft und zwar die „Größe“ zu simulieren (z.B. FPM, PARSIVAL).

Das Design der relevanten Produkteigenschaften von Nanopartikeln ist unmittelbar auf die avisierten Anwendungen zugeschnitten. Die Partikelgröße ist eine wichtige Eigenschaft, aber auch die chemische Zusammensetzung, die Kristallmodifikation, die Größe der kristallinen Bereiche innerhalb der Nanopartikel, die Anzahl der Primärpar-tikel innerhalb eines agglomerierten Partikels, die Primärpartikelgröße und auch der Sintergrad der Primärpartikel können eine wichtige Rolle in der Anwendung spielen. Besonders im Bereich der Funktionsanwendungen, wo die Nanopartikel aktiv wirken (z. B. als Katalysator oder als Funktionsschicht in der Gassensorik wobei die elektri-schen Eigenschaften eine Funktion der Gaszusammensetzung sind) sind Kenntnisse über diese Eigenschaften essentiell. Sehr oft werden Produkte empirisch entwickelt, da die Zusammenhänge zwischen den Produkteigenschaften der Nanopartikel und den Produkteigenschaften des Endproduktes nicht bekannt sind. Als Beispiel kann hier Industrieruß genannt werden, der als Füllstoff in der Gummiindustrie hauptsäch-lich für Autoreifen und Förderbänder verwendet wird, Hier ist die räumliche Struktur der Ruß-Agglomerate entscheidend, aber aufgrund fehlenden Grundverständnisses werden die Füllstoffe immer noch weitestgehend empirisch entwickelt. Ein weiteres Beispiel ist TiO2, eingesetzt als Pigment, dessen Farbe nicht nur eine Funktion der chemischen Zusammensetzung, sondern auch von der Partikelgröße und Partikelform ist.

Die Komplexität nimmt noch zu, da die genannten Eigenschaften oft eine Verteilungs-funktion aufweisen. Dies ist bekannt für die Partikelgröße, jedoch weniger für die Form, chemische Zusammensetzung, Ladungszustand, Sintergrad, Kristallmodifikati-on, etc. Die Partkelgröße wird bisher üblicherweise mittels Populationsbilanzen simu-liert, was aber für Systeme mit mehrdimensionalen Eigenschaften sehr zeitaufwändig ist und ein hohes Maß an Erfahrung voraussetzt. Ein an der Universität Duisburg-Essen entwickeltes Monte-Carlo-Simulationsverfahren erlaubt eine relativ einfache und numerisch effiziente Simulation verteilter Eigenschaften. Zwischen 2006 und


2008 wurde im Rahmen des AiF/IGF-Projekts "Einfluss von Ladungseffekten bei Na-no-Partikeln in Hochtemperaturverfahren" eine kombinierte Methode zur gekoppelten Simulation gasgetragener Partikelphasen erstmalig implementiert. Die Berechnung der Gasphase erfolgte hierbei noch grundsätzlich unter Verwendung des kommerziel-len Strömungslösers Fluent®. Für die Partikelphase wurde Monte-Carlo Verfahren an die spezifischen Erfordernisse angepasst. Eine einseitige Kopplung (sequentielle Be-rechnung von Kontinuumsströmung und Partikelphase) ergänzte die Arbeiten. Die Darstellung einer Partikelgrößenverteilung auf der Basis von z. Z. etwa 1000 Simulati-onspartikeln pro Zelle erlaubt die Berechnung von weiteren produktrelevanten Eigen-schaften.

1.2 Anlass für den Forschungsantrag

Ziel des Forschungsvorhabens ist die Bereitstellung von Methoden zur Simulation re-präsentativer Partikelpopulationen in gasgetragener Entstehung, vorwiegend in Pro-duktionsanlagen (Synthesereaktoren), mittelfristig jedoch auch für umweltrelevante Fragestellungen wie die Russ- und Staubbildung in natürlichen und technischen Pro-zessen. Der Anlass für den Forschungsantrag ergibt sich aus der Tatsache, dass der-artige Methoden erst ansatzweise und mit stark vereinfachter Modellbildung für indust-rierelevante Fragestellungen verfügbar sind. Da ferner ein experimenteller Zugang in vielen technischen Anlagen nur unter extremem Aufwand oder gar nicht möglich ist und eine neue Anlage i. A. erst nach Festlegung wesentlicher Konstruktionsmerkmale experimentell untersucht werden kann, sind entsprechende Simulationsmethoden un-erlässlich für Synthesereaktor- und Produktdesign.

Die Entwicklung und Anwendung der Software zur Auslegung und Optimierung von Reaktoren zur Herstellung von dispersen Produktpartikeln und Nanomaterialien wird erschwert durch die hohen Lizenzkosten kommerziell verfügbarer CFD-Software. Vor allem für kmU stellt dies eine in vielen Fällen unüberwindbare Hürde bei der Einfüh-rung von Simulationsverfahren in den Entwicklungsprozess dar.

Darüber hinaus wird die Integration zwischen der CFD- und der Monte-Carlo Methode, z. B. durch eine Zwei-Weg Kopplung, wobei die disperse Phase eine Rückkopplung auf die kontinuierliche Phase ausübt, auch nur möglich sein, wenn die CFD Software offen liegt, d.h. manipulierbar ist. Die äußerst anspruchsvolle Modellierung der Parti-kelbildungsprozesse erfordert in der Regel tiefe Eingriffe in die Programmstruktur be-stehender Software. Dies ist umso mehr von Bedeutung, da mit dem Vorhaben eine experimentell gestützte, sukzessive Modell- und Verfahrensvalidierung der Partikel-synthese angestrebt wird. Derzeit ist man bei der Planung und Auslegung von Syn-theseanlagen im Technikums- oder Pilotmaßstab oftmals nur auf Erfahrungswerte der Betreiber angewiesen.

Da die Monte-Carlo Software auf Linux-Clustern mit open source Software aufgebaut wurde, ist es nahe liegend, auch bei der CFD-Methode auf lizenzfreie Software zu-rückzugreifen. Eine Kopplung mittels Fluent® wäre auch aufgrund der erforderlichen großen Anzahl von PCs in dem LINUX-Cluster extrem kostspielig, da für jeden PC


Lizenzgebühren bezahlt werden müssen. Aus derzeit verfügbaren CFD-Programmen sticht das OpenFOAM Projekt durch eine hohe Qualität, gute Bedienbarkeit und eine große Anwendergemeinde hervor und bietet sich als quelloffenes (open source) Pro-gramm an.

2 Forschungsziel

Ziel des Forschungsvorhabens ist, wie bereits dargelegt, die Entwicklung von Pro-dukteigenschaften von nanopartikulären dispersen Systemen voraussagen zu können. Auf Grundlage des AiF/IGF-Projekts 14720 N wird die erfolgreich entwickelte ortsauf-gelöste Monte Carlo Simulationsmethode verwendet und weiter entwickelt. Die nume-rische Methode wird mittels ortsaufgelöster in situ Messungen der Produkteigen-schaftsverteilung von Nanopartikeln validiert. Es soll gezielt danach geforscht werden, wie Produkteigenschaften von Nanopartikeln in industriellen Prozessen entstehen. Mit diesem Verständnis wird es anhand von Simulationsergebnissen möglich sein, die jeweiligen Prozesse zu optimieren und wirtschaftlicher zu gestalten.

Ermöglicht werden soll dadurch:

• eine zuverlässigere nachträgliche Optimierung von Anlagen, in denen Produkt-partikel hergestellt werden,

• die Entwicklung neuer Produktionsmethoden basierend auf den Produkteigen-schaften,

und

• die Optimierung der Handhabungsproblematik, wie z. B. das Vermeiden von Partikelanlagerungen oder Ausnutzen der Partikelladung zur Partikelabschei-dung, durch eine gezielte Optimierung der Prozessbedingungen.

Die Qualität der entwickelten Methoden und Modelle wird anhand des am Institut für Verbrennung und Gasdynamik verwendeten wandbeheizten Rohrreaktors experimen-tell validiert. Insbesondere sollen hier einzigartige Diagnosemethoden zum Einsatz kommen, mit dem Ziel, den Partikelbildungsprozess in enger Wechselwirkung mit der Modellentwicklung besser verstehen zu können. Ohne dieses Verständnis ist das Upscaling im Sinne von Produktionsreaktoren nicht möglich.

3 Zusammenfassung der Ergebnisse

Kenntnisse über effiziente Methoden zur Einstellung der gewünschten Produkteigen-schaften in der Nanopartikelproduktion bedeuten für Firmen einen erheblichen Wett-bewerbsvorteil. Dabei werden Probleme mit der Partikelgrößen- und Produkteigen-schaftsverteilung oft durch räumliche Inhomogenitäten in der Strömung verursacht. Dadurch können zahlreiche Anwendungen für hochdefinierte Nanopartikel nicht reali-siert werden, was eine Verlangsamung der technologischen Entwicklung zur Folge hat.


Diese Probleme können mittels der in diesem Projekt realisierten Simulationsmethode analysiert werden. Die Anforderungen an die zu entwickelnde Methode entstanden zum einen aus den im Vorgängerprojekt erkannten Defiziten der bisherigen Software und zum anderen aus erhöhten Ansprüchen im Bezug auf Auflösung und Rechenzeit. So war die im Vorgängerprojekt benutzte Methodik für den Partikeltransport zu re-chenintensiv und hat zu Simulationszeiten von über eine Woche geführt. Darüber hin-aus war sie nur geeignet für fest vorgegeben Zellformen die mit den modernen un-strukturierten CFD-Gittern nicht kompatibel sind. Der traditionelle Programmierstil aus der Ingenieurspraxis hatte sich als unpraktisch herausgestellt und hatte zu der Forde-rung nach einem modularen Programmaufbau geführt. Vorteile sind eine leichte Adap-tierbarkeit auch für eine willkürliche Zahl von Partikeleigenschaften.

Da auf Basis von modernen und sehr preiswerten Grafikkarten mit vielen Hunderten Prozessoren mittlerweile auch anspruchsvolle numerische Simulationen durchgeführt werden können, wurde entschieden, die gesamte Software auf diesen sogenannten GPUs (graphics processing units) anzuwenden. Ein solches System ist auch für kmUs erschwinglich, die so auch mit wenigen PCs ein leistungsfähiges Rechensystem auf-bauen können. Die Verwendung von GPUs erforderte aber eine komplette Überarbei-tung der Algorithmen für die Partikeldynamik. Es wurde festgestellt, dass die bis dahin benutzte inverse Methode für ein System, wo zu jedem Zeitschritt viele neue Simulati-onspartikel in den Zellen eintreffen, ungeeignet ist. Dagegen zeigte sich, dass sich die A/R (acceptance/rejection) Methode bei großer Anzahl an Zellen auf der GPU effizient beschleunigen lässt, etwa um den Faktor 100. Da die Simulationspartikel auf der GPU Karte gespeichert sind und eine Datenübertragung zur CPU einige Minuten pro Zeit-schritt kostet, wurde auch die Berechnung des Partikeltransports mittels CUDA auf der GPU durchgeführt.

Neben der Berechnung von verteilten Partikeleigenschaften wurde im Rahmen des Vorhabens ein OpenFOAM Applicationsolver für die Aufgabenstellung angepasst, so dass die für die Monte Carlo Simulation notwendigen ortsaufgelösten Eingangsgrößen mit Hilfe von entsprechenden CFD-Simulationen ermittelt werden konnten. Dazu ge-hörten neben dem Geschwindigkeits- und Temperaturfeld auch Konzentrationsfelder sowie die Partikelbildung aufgrund homogener Nukleation. Zur numerischen Erfas-sung des chemischen Zerfalls des Prekursors wurde eine vereinfachte Reaktionskine-tik in OpenFOAM integriert. Aufgrund der Ortsauflösung der Größen in CFD-Verfahren ist die Erfassung komplexer Reaktionsmechanismen nur bedingt durchführbar und in jedem Fall mit einem signifikanten Anstieg in der Rechenzeit verbunden. Aufgrund der Erfahrungen aus dem Vorgängerprojekt hat es sich als sinnvoll herausgestellt, auch die Partikelbildung aufgrund homogener Nukleation im Rahmen der CFD-Simulation zu berechnen. Für die zur Beschreibung dieses Vorgangs notwendigen Gleichungs-systeme wurde eine OpenFOAM-Programmbibliothek entwickelt, die gekoppelt an einen Applicationsolver die Berechnung der zeit- und ortsaufgelösten Nukleationsrate ermöglicht. Um die durch die Partikelbildung auftretende Verarmung der Gasphase an der partikelbildenden Spezies zu erfassen, wurde ein Senkenterm in die Erhaltungs-


gleichung der Konzentration der partikelbildenden Spezies implementiert. Im Laufe des Projekts hat es sich als numerisch sinnvoll und effizient herausgestellt, die im Rahmen der Experimente gemessenen Wandtemperaturprofile als Temperaturrand-bedingung vorzugeben. Zu diesem Zweck wurde ein in OpenFOAM vorhandenes Mo-dul zur Vorgabe konstanter Temperaturwerte in Abhängigkeit einer Ortskoordinate adaptiert und um eine ebenfalls vorhandene Interpolationsroutine erweitert.

Es wurden Vergleichsrechnungen mit der kommerziellen CFD-Software Fluent durch-geführt, um den modifizierten OpenFOAM Applicationsolver zu validieren. Die Ergeb-nisse dieser Vergleichsrechnungen ergaben eine gute Übereinstimmung der Ergeb-nisse der beiden CFD-Verfahren und bestätigten die gute Eignung von OpenFOAM für die zu untersuchende Fragestellung.

Die im Rahmen der CFD-Simulationen ermittelten Daten für Strömungs- und Tempe-raturfelder sowie für die Konzentrationsfelder und die Partikelbildung durch homogene Nukleation wurden an das Monte-Carlo-Simulationsverfahren zur Berechnung der wei-teren partikeldynamischen Vorgänge übergeben. Die Methoden wurden validiert und auf ihre Eignung für den industriellen Einsatz getestet.

Hierzu wurden im Rahmen des Vorhabens eigene experimentelle Untersuchungen zur Nanopartikelsynthese durchgeführt. Anders als ursprünglich geplant, wurde hierbei kein durch Mikrowellen (MW) beheizter Reaktor verwendet. Die Erfahrung in anderen Projekten hat gezeigt, dass dieser Reaktortyp besonders schwer zu modellieren ist. Die Strömung ist nicht mit ausreichender Zuverlässigkeit zu bestimmen und somit auch der Wärmeeintrag, bzw. die Temperaturverteilung im Reaktor. Gleiches gilt für die Gasphasen-Kinetik des Prekursors. Obwohl also MW-beheizte Reaktoren imstan-de sind Nanopartikel in guter Qualität (bezüglich der experimentellen Reproduzierbar-keit und der Größenverteilung) herzustellen, sind die Unwägbarkeiten bezüglich der Modellierung zu groß um als Basis eines Validierungsexperimentes zu dienen.

Daher wurde als experimentelles Vergleichssystem die Kombination aus Heißwand-Reaktor, SMPS (scanning mobility particle sizer) und TEM (Transmissions-elektronenmikroskopie) als beste, realisierbare Variante ausgewählt. Als Stoffsystem wurde Eisen verwendet, das durch thermische Pyrolyse aus Eisenpentakarbonyl (Fe(CO)5) erzeugt wurde, da hierfür die Sicherheitsanforderungen moderat sind und langjährige Erfahrung mit diesem Prekursor vorliegt [1,2,3,4,5]. Als besondere Vorteile sind zudem der einfache Aufbau und die einfache Strömungsführung zu nennen, die zusammen mit den gut bekannten Materialeigenschaften der verwendeten Gase gu-ten Zugang zur CFD-Modellierung des Prozesses erlaubte.

Das Ziel des Vorhabens wurde erreicht.


4 Lösungsweg

4.1 Modellentwicklung

In bisherigen Forschungsarbeiten der Antragsteller wurde die Keimbildung und das Partikelwachstum mittels kommerzieller Software berechnet (Fluent®, FPM) und dann an die Monte-Carlo Software übergeben. Allerdings hat sich herausgestellt, dass das erforderliche Zusatzmodul zur Kalkulation der Konzentration des Precursormaterials, bspw. bei der Rußproduktion, nicht mit dem Fine Particle Model (FPM) kompatibel ist.

In der künftigen Anwendung der kombinierten CFD/Monte Carlo Software soll die Ab-hängigkeit von kommerzieller Software vermieden werden. Der Verbrauch an Precur-sormaterial aufgrund von Keimbildung wird daher mit einem neu zu entwickelnden Modul in OpenFOAM berechnet. Die neue integrierte Software ist im Erfolgsfall in der Lage, abwechselnd einen CFD Zeitschritt für die kontinuierliche Phase und einen Monte Carlo Zeitschritt für die disperse Phase auszugeführen.

Die bisher von FS 1 entwickelte Monte Carlo Software wurde spezifisch an wenige ausgewählte Fallbeispiele angepasst. Die neue Software ist durch einen modularen Aufbau allgemeiner zu gestalten. Hierzu soll diese aus verschiedenen Modulen mit definierten Schnittstellen aufgebaut werden. Da dies ein großes Maß an Erfahrung in Softwareentwicklung als auch Erfahrung mit CFD voraussetzt, wurde bei der Architek-tur dieser Software die Hilfe eines spezialisierten Software-Unternehmens (NUMRAX aus Duisburg) in Anspruch genommen. Die Software soll die verteilten Partikeleigen-schaften 1 bis N simulieren

4.1.1 Grundlegende Problemstellung

Ziel ist es, partikelbeladene Strömungen zu simulieren, die sich insbesondere durch starke Wechselwirkungen zwischen den Partikeln auszeichnen und hierbei ein breites physikalisches Spektrum abdecken. Aufgrund mehrdimensionaler Verteilung der Par-tikeleigenschaften, z.B. Partikelmasse, Partikeloberfläche, stoffliche Zusammenset-zung, oder weiterer Eigenschaften, scheidet eine Klasseneinteilung im Sinne eines sektionalen Modells aus. Grund hierfür ist, dass sich die je Eigenschaft erforderlichen Anzahlen der Sektionen multiplizieren. Wären beispielsweise nur 20 Sektionen für eine angemessene Beschreibung der Masse wie auch der Oberfläche erforderlich, so ergäbe sich bereits eine zweidimensionale Aufteilung in 20x20=400 Sektionen, die in einem Eulerschen Ansatz je Stützstelle eines Rechennetzes zu speichern und zu transportieren wären, wobei zusätzlich jede Sektion mit jeder anderen wechselwirkt. Es versteht sich, dass die numerische Komplexität, d. h. die erforderlichen numeri-schen Ressourcen als Funktion der Rechengenauigkeit bereits bei nur zwei Eigen-schaften eine hohe Potenz aufweist. Die Hinzunahme weiterer Eigenschaften ließe schließlich den Ressourcenbedarf förmlich explodieren.


Diese Einsicht führt zum Konzept der statistisch motivierten Monte-Carlo Simulation. Statt der Beschreibung der Mannigfaltigkeiten in einem fest vorgegebenen Variablen-raum, werden Simulationspartikel definiert, die jeweils die erforderliche Anzahl an Ei-genschaften tragen. Da reale Strömungen i. A. eine sehr große Partikelanzahldichte aufweisen, scheidet eine individuelle Berechnung einzelner Partikel aus. Vielmehr sol-len die Simulationspartikel je für ein Ensemble an realen Partikeln stehen, so dass die numerisch zu prozessierende Anzahl an Simulationspartikeln beherrschbar bleibt.

Die physikalischen Vorgänge in partikelbeladenen Gasströmungen lassen sich in fol-gende Teilprozesse einordnen:

• Partikeltransport

• Partikelquellen und ggf. Senken

• Wechselwirkung der Partikel untereinander

▪ Optional zusätzlich: Wechselwirkung der Partikel mit der umgebenden Gasphase

Die numerische Modellierung orientiert sich grundsätzlich auch an dieser Einteilung. Eine objektorientierte Programmierung ermöglicht hierbei die modulare Anwahl ver-schiedener physikalischer Detaillierungen und numerischer Ansätze.

Obige Reihenfolge ist zunächst irreführend, da ein Partikel erst entstehen muss, bevor es transportiert werden kann. Die Modellierung der Transportterme ist jedoch für die Implementierung der übrigen numerischen Modelle von großer Bedeutung, weswegen der Partikeltransport und somit das Transportmodell vorgezogen werden.

4.1.1.1 Präambel zum Rechenleistungsbedarf und numerischen Aufwand

Voran stehend wurde bereits dargestellt, dass eine klassische, sektionale Methode zur Simulation von Partikelpopulationen mit vielfältigen Eigenschaften sehr schnell an ihre Grenzen stößt, was zugleich die wesentliche Motivation für die gewählte Monte-Carlo Methode ist. Obwohl letztere eine im Vergleich zu einem sektionalen Ansatz deutlich effektivere numerische Umsetzung und somit eine wesentlich engere Begren-zung des Variablensatzes ermöglicht, handelt es sich weiterhin um außergewöhnlich mächtige Problemstellungen. Auch bei Annahme zukünftig weiter sinkender Kosten für die erforderliche Rechenleistung, kann aber grundsätzlich noch nicht von der voll-ständigen Deckung des Bedarfs an Rechenleistung gesprochen werden. Vielmehr ist die Aussagekraft jeder Simulation durch die verfügbare Rechenleistung begrenzt. Aufgrund einer relativ niedrigen numerischen Komplexität der Monte-Carlo Simulation, werden die prognostizierten Leistungssteigerungen zukünftiger Rechenanlagen die dann mögliche Aussagekraft der Simulationen jedoch erheblich steigern.


Zur Erläuterung des hohen Leistungsbedarfs soll der Vergleich zu klassischen, strö-mungsmechanischen Simulationen angestellt werden: So erfordert eine 3D-Navier-Stokes Simulation inklusive Turbulenzmodellierung typischerweise sechs bzw. sieben Variablen je Stützstelle eines Gitters für eine inkompressible, respektive kompressible Strömung. Auch unter Hinzunahme chemischer Spezies, z. B. für die Simulation von Flammen oder anderweitig reagierender Strömungen, bleibt der Variablensatz je Stützstelle des Gitters auf ca. ein bis zwei Dutzend begrenzt. Um die Mannigfaltigkei-ten der Eigenschaften einer Partikelpopulation hinreichend präzise beschreiben zu können, sind im Gegenzug tausende Simulationspartikel je Stützstelle eines Gitters erforderlich.

Ähnliche Extreme gelten für die räumliche und zeitliche Diskretisierung. Die Orts-diskretisierung einer einfachen Strömungssimulation erfolgt üblicherweise mittels ei-nes Eulerschen Ansatzes für die Transportterme. Die Zeitintegration schließlich erfolgt über explizite bzw. besser implizite Ansätze und erlaubt Zeitschritte (d. h. zeitliche Fortschritte der Simulation), die eine Signalausbreitung in der räumlichen Größenord-nung einer Gitterzelle zulassen. Bei Partikelpopulationen sind demgegenüber je Zeit-schritt zunächst alle denkbaren Wechselwirkungen aller Partikel zu berücksichtigen, die sich jedoch über stochastische Ansätze (namensgebend für die Monte-Carlo Si-mulation) begrenzen lassen. Sollen diffusive Effekte ausreichend berücksichtigt wer-den, sind aufgrund der hohen Mobilität kleiner Partikel zudem nur extrem kleine Zeit-schritte zulässig.

Zusammenfassend erfordert die Simulation partikelbeladener Strömungen im Ver-gleich zu einer klassischen Strömungssimulation

• deutlich größere Variablensätze,

• bei gleichzeitig wesentlich stärkerer Interaktion

• und ebenso wesentlich kleinere, zulässige Zeitschritte.

Obige Eigenschaften der detaillierten Simulation einer suspendierten Partikelpopulati-on begrenzten aufgrund des hohen Rechenleistungsbedarfs in der Vergangenheit die Anwendung auf rein akademische Testfälle. So wurden zunächst 0-dimensional-instationäre Modelle zur Simulation des zeitlichen Verlaufs einer Partikelpopulation entwickelt. Um eine räumliche Aussage zu ermöglichen, konnten die Ergebnisse die-ser 0D-Simulationen längs der Trajektorien entlang der Stromlinien, z.B. in einem Par-tikelreaktor, aufgetragen werden. In Sonderfällen ermöglichte dies eine ein-, zwei- oder gar dreidimensionale stationäre, d. h. nicht mehr zeitabhängige Simulation. Ne-ben der Beschränkung auf stationäre Simulationen, konnten mit dieser Methode Prob-lemstellungen, bei denen eine Wechselwirkungen mit der Gasphase, eine ausgepräg-te diffusive Durchmischung der Partikel selbst und/oder Rückströmgebiete eine we-sentliche Rolle spielen, kaum erschlossen werden.


4.1.1.2 Begründung für die Benützung von Grafikkarten für die Simulation von

praxisnahen Reaktoren

Industrielle Aerosol-Reaktoren sind konstruktiv gesehen sehr komplexe dreidimensio-nale (3D) Systeme, auf die oben genannte Einschränkungen nicht zutreffen und die somit nur in seltenen Fällen mittels 0D trajektorienbasierten Lösungen zu beschreiben sind. Realitätsnahe Simulationen erfordern daher Rechengitter, die ähnliche Anzahlen an Stützstellen (Gitterpunkte oder Zellen je nach Verfahren) aufweisen, wie dies bei der Simulation reiner Fluidströmungen der Fall ist. Idealerweise sollten die gleichen Rechengitter für die Simulation der Gasströmung und der ortsaufgelösten Partikelpo-pulationen Verwendung finden. Die erforderliche Zellenzahl liegt dann in der Größen-ordnung von ca. 1 Mio. Stützstellen, exemplarisch etwa für ein Gitter aus 100x100x100 Stützstellen. Zum Vergleich hierzu sei erwähnt, dass das im Vorläufer-projekt (AiF/IGF 14720N) entwickelte Verfahren, ebenfalls auf Basis einer ortsaufge-lösten Monte-Carlo Simulation, neben weiteren Einschränkungen bereits für nur ca. 1000 Zellen bei knapper Anzahl an Simulationspartikeln Rechenzeiten im Bereich ei-niger Tage erforderte.

Die vollständig dreidimensionale, instationäre und detaillierte Simulation einer parti-kelbeladenen Strömung stellt somit eine außerordentlich herausfordernde Aufgabe dar. Es versteht sich daher, dass nur die Kombination aus maximaler numerischer Effizienz einem hohen Grad an Parallelisierung und maximaler Wirtschaftlichkeit der notwendigen Rechnerinfrastruktur Erfolg verspricht.

Es sei an dieser Stelle vorweggenommen, dass die Simulation für Parallelrechner mit abgesetzten Stream-Computing Einheiten unter CUDA optimiert wurde. Hierbei han-delt es sich technisch um Grafikkarten (GPUs: graphics procesing units), die intern je Einzelstück bereits aus massiv parallelen Prozessierungseinheiten bestehen. Diese Technik ermöglicht die derzeit mit Abstand günstigsten Preis-Leistungsverhältnisse und ebenfalls sehr günstige Stromverbrauchswerte. Letztere sind wirtschaftlich bedeu-tend: Werden preiswerte Komponenten des Massenmarkts (PC-Processoren, Grafik-karten, Speichersysteme) verwendet und zu parallelen Großrechenanlagen zusam-mengefügt, übersteigen die auf die Nutzungsdauer hochgerechneten Betriebskosten bereits die Investitionskosten, Tendenz weiter steigend. Den wirtschaftlich günstigen Bedingungen GPU-basierter Rechencluster steht allerdings auch eine sehr diffizile und empfindliche Programmoptimierung gegenüber, die nicht nur detaillierte Algorith-menteile, sondern auch bereits das grundlegende Datenhaltungskonzept der vorlie-genden Entwicklungen maßgebend prägt.

Erschwert wird dieser informatisch-algorithmische Teil der Entwicklungen zudem durch die geometrischen Anforderungen an eine industriell einsetzbare Simulations-software. Nahezu alle kommerziell oder frei verfügbaren Strömungslöser basieren auf unstrukturierten oder hybriden Gittern oder sie besitzen zumindest die grundsätzlichen Fähigkeiten auf diesen Gittertypen lauffähig zu sein. Dies aus gutem Grund: Techni-


sche Apparate weisen häufig komplexe Geometrien auf, die durch strukturierte, soge-nannt IJK-abzählbare Gitter nur schwer vernetzt werden können. Selbstverständlich wird daher auch für die Monte-Carlo Simulation die Lauffähigkeit auf unstrukturierten und hybriden Gittern verlangt.

Obige Anforderung maximaler geometrischer Flexibilität steht in einem Zielkonflikt zu einer Optimierung für GPU-basierte Rechnersysteme. Grund hierfür ist, dass die mas-siv parallelen GPU-Rechenwerke sehr empfindlich auf Adresssprünge reagieren und ihre überragende Leistungsfähigkeit nur ausspielen können, falls eine hohe Anzahl unabhängiger Rechenoperationen bei gleichzeitig hoher Rechenintensität gefordert ist. Wesentliche Teile der nachfolgend beschriebenen Entwicklung befassen sich mit dieser Problematik.

4.1.1.3 Vorgehensweise

Die wichtigsten Anforderungen an die Software, und wesentliche Änderungen gegen-über der alten Simulationsmethode, sind folgende.

1. Benutzung von Gittern beliebiger Struktur

Heutzutage ist dies der Standard. Dies ist wichtig für die Kompatibilität zur open sour-ce CFD Software wie OpenFOAM.

2. Beschreibung willkürlicher Partikeleigenschaften

Die Software aus dem Vorgängerprojekt war nicht modular aufgebaut und erlaubte keine einfache Erweiterung der Anzahl der Partikeleigenschaften. Der Grund dafür war die komplexe Datenübertragung zwischen den einzelnen Simulationsdomänen auf den Clusterknoten mittels MPI. Eine komplette Programmänderung wäre notwen-dig gewesen, dies war aufgrund der Größe (etwa 150 Seiten Quellcode) zu aufwen-dig.

3. Verwendung von open-source Software

Die Strömungssimulationen (CFD) sollen mit einer open-source Software durchgeführt werden. Wesentlich ist hier die Möglichkeit, die Software genau kennenzulernen und anzupassen da der Quellcode offen ist. Diese Anforderung ist maßgeblich für den Einsatz der Software bei KMUs. Für das Vorhaben wurde die kostenlose CFD Soft-ware OpenFOAM benutzt. OpenFOAM ist in allen drei Forschungsstellen im Einsatz und hat sich in den letzten Jahren als Standard für open source CFD Software etab-liert.

4. Vereinfachung des Transportes von Simulationspartikeln


Im Vorgängerprojekt basierte der Partikeltransport auf der Wahrscheinlichkeit, dass ein Partikel innerhalb eines Zeitschritts t von der Position 1 in der Zelle 1 in die Positi-on 2 in der Zelle 2 transportiert wird. Diese Wahrscheinlichkeit ist z. B. in Zylinderko-ordinaten:

(((( ))))

−−−−−−−−++++−−−−====

−−−−====

Dt

rrrr

DtDt

PP

DtPd

4

cos2exp

4

1

4exp

4

1 21212

22

1

2

21 θθ

ππ (1)

wobei D der partikelgrößenabhängige Diffusionskoeffizient ist. Die relevanten geomet-rischen Parameter sind in der Abbildung 4.1 dargestellt.

•

Abbildung 4.1: Partikeltransport in Zylinderkoordinaten.

Die Wahrscheinlichkeit dass ein Partikel sich von der Zelle 1 in die Zelle 2 bewegt wird dann mittels numerischer Lösung folgenden Integrals bestimmt, d. h. über alle mögli-chen Kombinationen der Positionen 1 und 2.

212121 θθ dddrdrrrPP d∫∫∫∫ ∫∫∫∫ ∫∫∫∫ ∫∫∫∫==== (2) Der Zeitschritt wird so ausgewählt, dass der Transport in die übernächsten Nachbar-zellen minimal ist. Durch die Verwendung von variablen Gewichtungsfaktoren, die die Partikelkonzentration darstellen, wird in allen Zellen eine ähnliche Statistik gewährleis-tet. Diese Methode ist präzise, aber erstens sehr zeitintensiv (Vielfachintegral !): die Rechenzeit lag bei etwa einer Woche für 1000 Zellen x 1000 Simulationspartikel, An-zahl der Zeitschritte war etwa 1000. Zweitens ist sie nur praktikabel bei fest vorgege-ben Zellformen. Die Benutzung von unstrukturierten Gittern ist damit nicht möglich.

5. Beschleunigung der Simulation partikeldynamischer Vorgänge mittels der Monte Carlo Technik.

Die Programmstruktur erfordert eine Abwechslung zwischen Transportschritt und Par-tikeldynamikschritt. Es gibt viele Zeitschritte, um zu verhindern, dass die Simulations-partikel eine Zelle überspringen. Darüber hinaus muss dann die Partikeldynamik in vielen tausend Zellen simuliert werden. Bei 106 Zellen und 104 Zeitschritten würde dies bedeuten, dass 1010 mal ein partikeldynamischer Simulationsschritt durchgeführt


werden muss. Das heißt, dass die Monte Carlo Technik für die Partikeldynamik stark beschleunigt werden muss, da im Vorgängerprojekt einige Sekunden pro Koagulati-ons-Zeitschritt und pro Zelle benötigt wurden – dies würde eine Simulationszeit von etwa 1000 Jahren erfordern.

Diese fünf wesentlichen Anforderungen waren die Motivation der durchgeführten Ar-beiten.

4.1.2 Optimierung der Monte Carlo Simulationstechnik für die Partikeldynamik

Die ersten Bemühungen der Forschungsstelle 1 die Monte Carlo Simulation zu be-schleunigen basierten auf der inversen Methode der Monte Carlo Technik (Kruis et al., 2000). Diese Methode ist schematisch in Abbildung 4.2 dargestellt.

Abbildung 4.2: Inverse Methode für die Auswahl eines Partikelpaares in der Koagulation.

In dieser Methode wird ein Partikelpaar (i,j) ausgewählt, wenn

′Ckk=1

i−1

∑ ∆t / 2 < r < ′Ck∆t / 2k=1

i

∑ , i ∈ 1, N[ ]

(3)

und

′βkk=1

j−1

∑ ∆t / 2 < ∆r < ′βk∆t / 2k=1

j

∑ , j ∈ 1, N[ ]

(4)

wobei r eine Zufallszahl zwischen 0 und 1 ist, und


∆r = r − ′Ckk=1

i−1

∑ ∆t 2 (5)

′Ci ist definiert als Gesamtkoagulationsrate ′Ci

′Ci =1

V2′βij

j=1, j≠i

N

∑

(6)

wobei V das Systemvolumen ist, N die Gesamtzahl der Simulationspartikel, ′βij der

modifizierte Koagulationskoeffizient zwischen Partikel i und j. Diese Methode lässt sich sehr gut und effizient auf der GPU implementieren, da hier sehr viele einfache Summen berechnet werden müssen. Die erreichten Beschleunigungsfaktoren sind beträchtlich und hängen von der Zahl der Simulationspartikel ab (Abbildung 4.3). Aus dem Bild wird klar, dass der wichtigster Faktor die Anzahl der betrachteten Simulati-onspartikel ist: je größer die Zahl, umso größer die Beschleunigung. Hierbei ist zu be-achten, dass bei der benutzten Grafikkarte GTX285 die Anzahl der Prozessoren 256 beträgt. Es hat sich herausgestellt, dass die Benutzung einer GPU optimal ist, wenn einfache Operationen sehr häufig wiederholt werden.

Abbildung 4.3: Gemessene Beschleunigungsfaktoren (speedup) auf einer GPU (Modell GTX 285 mit 256 Prozessoren) im Vergleich zur CPU (1 Prozessor) als Funktion der dimensionslosen Koagulationszeit, für unterschiedliche Anzahl der Simulationspartikel.


Es hat sich aber auch herausgestellt, dass die Methode zwei wesentliche Probleme hat. Das erste Problem entsteht dadurch, dass die Methode nach jeder Kollision einen

Zeitschritt ∆t berechnet:

∆t =V

′Ckk=1

N

∑

Man hat in dieser Methode keine Kontrolle über den Zeitschritt, die Kollisionsfrequen-zen führen zu eindeutigen aber nicht beeinflussbaren Zeitschritten. Das zweite Prob-lem ist, dass diese Methode nur dann effizient ist, wenn sich das Partikelensemble nicht ändert. Dann kann nämlich nach jeder Kollision (Paarauswahl) eine einfache Korrektur der ′Ci Werte durchgeführt werden, und es müssen nicht alle Paarinterakti-

onen berechnet werden (‚smart bookkeeping‘). Das Problem entsteht, wenn durch Partikeltransport viele neue Partikel in einer Zelle hinzu kommen. Das macht es un-

möglich, das ‚smart bookkeeping‘ durchzuführen und es müssen alle 2N Kollisions-frequenzen und die Aufsummierungen berechnet werden. Aus dem Grund wurde dann auf die „Acceptance-rejection“ (A/R) Technik zurückgegriffen.

Die A/R Technik funktioniert folgendermaßen: Es werden zwei Simulationspartikel i

und j über Zufallszahlen selektiert. Die Kollisionswahrscheinlichkeit P wird berechnet nach

( )

max

,

β

β jiP =

(7)

wobei maxβ der maximale Koagulationskoeffizient der Population ist. Dann wird mittels

einer dritten Zufallszahl r dieses Partikelpaar akzeptiert (wenn Pr < ) oder abgelehnt (wenn Pr ≥ ). Danach wird der Zeitschritt bestimmt mittels:

mean

N

i

N

ijj

ij

N ββ

τ′

=

′

=

∑ ∑= ≠=

2

1 ,1

11

(8)

Auf dem ersten Blick scheint die A/R Methode für eine massiv parallele Berechnung ungeeignet zu sein, da man so lange Partikelpaare überprüft bis eine Überprüfung zu einem ‚accept‘ führt:

RRRRRRRRRRRRRRRRRRRRRRRRRA

wonach man aufhört und den Zeitschritt berechnet. Allerdings kann man auch die A/R Überprüfung massiv parallel machen lassen, und danach ein akzeptiertes Partikelpaar

auswählen. Dieser Ansatz wurde hier verfolgt. Ein weiteres Problem ist, dass maxβ

bekannt sein muss, um die Überprüfung durchzuführen. Frühere Arbeiten und andere

Arbeitsgruppen benutzen häufig einen Schätzwert für maxβ . So lange dieser Wert nicht

kleiner ist als der eigentlich korrekte Wert (basierend auf allen Partikelpaaren), führt


dies nicht zu einem Fehler, aber es führt zu einer zu großen „Reject“-Rate. Dies wird in Abbildung 4.4 Gezeigt: Mit Verlauf der Zeit nimmt diese Rate immer mehr zu, und die Simulation wird sehr langsam. Aus dem Grund haben wir untersucht, ob es mög-

lich ist, aus dem Schätzwert für die mittlere Koagulationsrate meanβ ′ , die aus den zufäl-

lig ausgewählten Partikelpaaren ermittelt werden kann, maxβ abzuschätzen. In

Abbildung 4.5 wird die zeitliche Entwicklung dieser beiden Koeffizienten für einen typi-schen Koagulationsverlauf gezeigt. Es zeigt sich ein etwa konstanter Wert γ bei

meanβγβ ′=max . Aufgrund der gleichbleibenden Form der Verteilungsfunktion (annähernd

selbsterhaltend) erscheint dies auch logisch.

Abbildung 4.4: Anzahl der nicht-akzeptierten Partikelpaar-Kombinationen als Funktion der dimensions-

losen Koagulationszeit, bei gleich bleibendem maxβ


Abbildung 4.5: Verlauf meanβ ′ und maxβ von als Funktion der dimensionslosen Koagulationszeit.

Allerdings ist nicht klar, ob man diesen Ansatz auch verfolgen darf, wenn die Vertei-lungsfunktion deutlich breiter oder sogar bimodal, wie in keimbildenden Systemen ist. Dies wurde anhand künstlich erzeugter Partikelgrößenverteilungsfunktionen unter-sucht. Abbildung 4.6 zeigt, dass die Konstante γ größer wird bei zunehmender Vertei-

lungsbreite, die Zunahme hält sich allerdings noch in Grenzen.


Abbildung 4.6: Konstante γ als Funktion der Verteilungsbreite (geometrische Standardabweichung).

Für Systeme wo homogene Keimbildung auftritt, entstehen bimodale Verteilungen. Dies wurde in Abbildung 4.7 simuliert. Auch hier ist der Wert von γ nicht allzu groß.


Abbildung 4.7: Konstante γ bei bimodalen Größenverteilungen, für verschiedene geometrische

Durchmesser und relative Partikelanzahlkonzentrationen.

Für die restliche Untersuchungen wurde ein Wert von γ = 40 angesetzt. Dieser ist

zwar größer als für selbsterhaltende Verteilungen notwendig, aber scheint die meisten hier betrachteten Fälle abzudecken und man baut damit eine gewisse Sicherheit ein. In Abbildung 4.8 ist die Struktur der parallelisierten Software dargestellt. In verschie-denen „Threads“ (kleine Programme die auf den Einzelprozessoren der GPU durchge-führt werden) werden Partikelpaare i ausgewählt und die zugehörigen Koagulati-

onskoeffizienten iβ berechnet. Dann wird die Summe ∑ iβ berechnet, dies erfordert bei 12 −n Threads nur n Rechenschritte, z. B 10 Schritte bei 512 Threads. Auf einer CPU hätte diese Summierung 512 Rechenschritte benötigt. Wenn alle Aufsummie-

rungen bereit sind (‚sync‘), wird das maxβ mit der bereits skizzierten Methode abge-

schätzt. Dies erlaubt dann für alle (z.B. 512) Partikelpaare die Berechnung der Kollisi-onswahrscheinlichkeit und mittels einer Zufallszahl die Entscheidung ob „A oder R“. Es wird dann ein „A“ Paar ausgewählt wenn alle Threads so weit sind (sync).


Abbildung 4.8: Struktur der parallelisierten Software, programmiert mittels CUDA.

Diese Methode wurde so implementiert und getestet. Es stellte sich heraus, dass die Beschleunigung durch die Verwendung der GPU nur dann eintritt, wenn das Pro-gramm massiv parallel betrieben wird, d. h. für viele verschiedene Partikelsammlun-gen. In unserem Fall können diese Partikelsammlungen Zellen darstellen. Dies ist in Abbildung 4.9 erkennbar. Es wurden Simulationen für 1,10, 100 und 1000 Zellen durchgeführt. Für N=1 ist die Verwendung von CUDA kaum begründbar, da sie für kleinere Simulationszeiten sogar langsamer ist als mit der CPU. Nur für größere Zel-lenzahlen gibt es eine erkennbare klare Beschleunigung, bis zu einem Faktor 20. Auch wird klar, dass die Erhöhung von 100 auf 1000 Zellen wenig Beschleunigung mehr bringt, dies ist angesichts der beschränkten Anzahl von GPU-Prozessoren (hier 256) auch verständlich.


Abbildung 4.9: Beschleunigungsfaktor durch den Einsatz einer GPU, für verschiedene Zellenzahlen (N) mit je 2000 Simulationspartikeln (GPU: NVIDIA GTX285 mit 256 Prozessoren)

Wichtige Schlussfolgerungen sind also:

1. Die Inverse Methode ist, durch die notwendige Initialisierung über alle Partikel-

paare, eine 2N Operation, für die gekoppelte Beschreibung der Partikeldyna-mik und Partikeltransport ungeeignet.

2. Nur die A/R Methode ist in der Lage, die nach jedem Transportschritt neue Par-tikelpopulation zu berücksichtigen.

3. Die A/R Methode ist nur bei großer Anzahl an Partikelsammlungen (hier: Zel-len) durch GPU-Parallelisierung bedeutend schneller als eine CPU. Allerdings ist das der hier interessierende Fall, da viele tausend Zellen gleichzeitig simu-liert werden sollen.

Die Software wurde noch weiter optimiert, z. B. durch Benutzung des schnelleren Speichers auf der GPU-Karte. Der CUDA-Programmierung ist komplex und hochspe-zifisch, so gibt es 7 unterschiedliche Speicharten auf der GPU-Karte, alle mit spezifi-schen Vor- und Nachteilen. Die benötigte Rechenzeit für eine unterschiedliche Anzahl an Simulationspartikeln pro Zelle, sowohl für den CPU als auch für die parallelisierte Software auf einer GPU, ist in Abbildung 4.10 ersichtlich. Die Beschleunigung auf-grund der GPU-Benutzung liegt bei etwa einem Faktor von 150.


Abbildung 4.10: Simulationszeiten für verschiedene Simulationspartikelzahlen pro Zelle für CPU und GPU, für 1000 Zellen.

4.1.3 Partikeltransportmodul für die Monte Carlo-Simulation auf Gittern beliebi-

ger Struktur

Die Partikelpopulation wird physikalisch durch vier Mechanismen transportiert:

• Transport aufgrund der Trägerströmung

• Stochastische, thermische Eigenbewegung

• Optional: Wirkung äußerer Kräfte, z. B. elektrische Kraftfelder

• Optional: Eigene Trägheit. I. A. nur für relativ große Partikel von Bedeutung

In vorliegender Entwicklung wurde zunächst nur der Transport aufgrund der Träger-strömung und die thermische Eigenbewegung (Diffusion) berücksichtigt. Dieses sind i. A. die wesentlichen Transportmechanismen, die in Partikelreaktoren eine Rolle spie-len.

4.1.3.1 Eulerscher versus Lagrangescher Partikeltransport

Die numerische Modellierung kann anhand zweier grundlegend unterschiedlicher An-sätze erfolgen, die hier gegenübergestellt werden sollen.


• Eulerscher Ansatz (Abbildung 4.11)

• Lagrangescher Ansatz (Abbildung 4.12)

Abbildung 4.11: Eulerscher Transport: Die Partikel können über die Zellwände hinweg in benachbarte

Zellen transportiert werden.

Abbildung 4.12: Lagrangescher Transport: Die Partikel werden entlang ihrer Bewegungstrajektorien verfolgt.

In beiden Ansätzen werden die Partikelpopulationen je Stützstelle des Netzes, im vor-liegenden Fall die Zellen, gehalten. Zu einer Zelle Z gehört somit eine Menge an Si-mulationspartikeln S(Z). Die Mächtigkeit von S(Z) liegt in einem Größenordnungsbe-reich von 103 bis 104 je Zelle.

4.1.3.2 Eulerscher Partikeltransport

Der Eulersche Transportansatz geht davon aus, dass jedes Simulationspartikel mit entsprechenden Wahrscheinlichkeiten entweder in der betreffenden Zelle verbleibt oder aber in eine benachbarte Zelle transportiert wird. Die exakte Position des Simula-


tionspartikels innerhalb der Zelle spielt hierbei keine Rolle. Wesentlich ist ausschließ-lich die Zugehörigkeit zu einer Zelle Z. Da die diffusive Mobilität der Partikel aufgrund der Überlagerung mit der Trägerströmung nicht aus bereits erwähntem 4-fach Integral herausgezogen werden kann, müssen die Wahrscheinlichkeiten je Partikel berechnet werden. Da mehrere Zellnachbarschaften berücksichtigt werden müssen, ist zudem die Bestimmung der Übergangswahrscheinlichkeiten für jede dieser Nachbarschaften erforderlich. Diesem hohen Rechenaufwand steht dann gerade einmal ein möglicher Partikelübertritt gegenüber.

In der vorliegenden Entwicklung sollen zudem 3-D Problemstellungen erfasst werden, d. h. die Integrale müssen über je drei Dimensionen für die möglichen Orte innerhalb des Start- und Zielgebiets geführt werden. Ferner sollen verzerrte, unstrukturierte und hybride Gitter zum Einsatz kommen, so dass einerseits die Zahl der möglichen Zell-nachbarschaften stark steigt und zudem die Start- und Zielgebiete geometrisch kom-plexe Gebilde werden, die sinnvoller Weise nur durch Tetraeder-Parzellierungen er-fasst werden könnten. Schließlich steigt die Anzahl der Zellen in 3-D Simulation ge-genüber 2-D Ansätzen sehr stark an und aus Gründen der Ergebnisgenauigkeit sollte im Vergleich zum Vorgängerprojekt mit deutlich größeren Partikelanzahlen je Zelle simuliert werden können.

Obwohl im Vorgängerprojekt gute Ergebnisse mit dem Eulerschen Verfahren erzielt werden konnten, ist diese Methode für die vorliegenden Entwicklungen zu rechenzeit-aufwendig. Es ließ sich abschätzen, dass bereits moderat aufgelöste Simulationen auf aktuell verfügbaren Rechenanlagen zu Rechenzeiten im Bereich Jahren bis Jahrhun-derten führen würden, ein völlig inakzeptabler Gedanke.

Die Eulersche Transportmethode wurde daher zugunsten einer Lagrangeschen Vari-ante verworfen.

4.1.3.3 Lagrangescher Partikeltransport

In Lagrangescher Formulierung trägt jedes Simulationspartikel räumliche Koordinaten, es ist also seine exakte Position bekannt. Die Zugehörigkeit zu einer Zelle des Netzes ergibt sich aus seiner Position. Da die Zellen den Raum eindeutig parzellieren, ist je-der Ort im Koordinatensystem und somit das Simulationspartikel auch eindeutig einer Zelle zuzuordnen. Die Bewegung der Simulationspartikel ergibt sich aus Superposition der beteiligten Transportmechanismen, wobei die Stochastik der Diffusion unmittelbar in die Berechnung individueller Trajektorien eingeht. Spielt die eigene Trägheit keine Rolle, entfällt die zeitliche Integration der Beschleunigung. Eine Integration der Ge-schwindigkeiten ist dann ausreichend, um die Wegstrecke und somit die Trajektorie eines Simulationspartikels festzulegen.

Bezüglich der Stochastik des diffusiven Transports, ist die Übertragbarkeit des indivi-duellen Verhaltens physikalischer Partikel auf das Verhalten der Simulationspartikel


nicht unmittelbar einsichtig. Diese Fragestellung stellte sich in gleicher Weise auch bei einer Eulerschen Betrachtung: Aufgrund der Zusammenfassung von Gruppen physi-kalischer Partikel zu Simulationspartikeln ergibt sich in einem einzelnen Zeitschritt ei-ne zu geringe Durchmischung, da das betreffende Simulationspartikel je nur einen repräsentativen Ort einnehmen kann. Da jedoch auch die weitere Erlebnishistorie der Partikel, vor allem die Koagulation mit anderen Partikeln, einer Stochastik unterworfen ist und da eine hohe Anzahl an Zeitschritten ausgeführt wird, verliert sich der Effekt entlang der zeitlichen Entwicklung. Die folgenden Abbildungen verdeutlichen dies schematisch: Abbildung 4.13 zeigt den Transport physikalischer Partikel. In Abbildung 4.14 sind diese als Simulationspartikel zusammengefasst. Abbildung 4.15 schließlich zeigt eine Menge an Simulationspartikeln längs ihrer konvektiv-diffusiven Trajektorien.

Abbildung 4.13: Mikroskopischer Transport physikalischer Partikel mittels Trägerströmung (determinis-tisch) und Diffusion (stochastisch).

Abbildung 4.14: Trajektorie eines diskreten Simulationspartikels als Ersatzmodell für die Partikelgrupp-pe.

Abbildung 4.15: Mannigfaltige Simulationspartikel mit stochastischer Bewegungsüberlagerung entlang ihrer Trajektorien.


Die Bewegung der Partikel ist sehr allgemein über deren Kinetik bestimmt. Die Bewe-gungsgleichungen der repräsentativen, diskreten Simulationspartikel lauten in allge-meiner Form:

∫

∫∆+=∆

∆+=∆

−=

Störung

Störung

sdtvs

vdtav

Fa

Hierbei sind a die Beschleunigung, F die auf das physikalische Partikel einwirkende

Kraft, v∆ die Geschwindigkeitsänderung und s∆ die Änderung des Ortsvektors längs

der Bewegungstrajektorie. Störungv∆ und Störungs∆ bezeichnen Störkomponenten. Da die

physikalischen Partikel in Gruppen zu Simulationspartikeln zusammengefasst werden, gelten die Gleichungen in analoger Form auch für diese, wobei die Beschleunigung weiterhin die eines einzelnen physikalischen Partikels darstellt.

Ausgehend von einem Anfangsort und einer Anfangsgeschwindigkeit kann somit der neue Ort des Partikels nach einem Zeitschritt bestimmt werden.

Die stochastische Brownsche Beschleunigung wird mittels einer Brownschen Zufalls-kraft modelliert:

( )t

SGtF ii

∆= 0

π

(9)

wobei iG eine Gauss-verteilte Zufallszahl ist und

cp

B

CSd

TS

2520

216

ρπ

νκ=

(10)

Obiger Ansatz geht von einer eigenen Kinetik der Partikel gegenüber der Trägerströ-mung aus. Eine Entkopplung der Diffusionsgeschwindigkeit ist erforderlich, jedoch ist die Anpassung an die konvektive Geschwindigkeit der Trägerströmung hinreichend schnell, um auf eine vollständig entkoppelte Kinetik verzichten zu können. Entfällt die eigene Kinetik, kann ebenfalls auf eigene Geschwindigkeitsfelder für die Simulations-partikel verzichtet werden, was aufgrund der hohen Partikelanzahl sowohl bzgl. der Rechenzeit als auch bzgl. der erforderlichen Speichermenge sehr vorteilhaft ist. Diese Vereinfachung ist in der gegenwärtigen Version des Verfahrens noch nicht verfügbar, wird jedoch zu erheblichen Leistungssteigerungen führen.


Es lässt sich eine Referenzgeschwindigkeit der Partikel definieren, die im statistischen Mittel des stochastischen Simulationsvorgangs zu einer äquivanten Diffusivität der Partikel führt.

Diese Referenzgeschwindigkeit hängt dann in guter Näherung nur noch von der Parti-kelmasse und den örtlichen Zustandsgrößen der Trägerströmung ab und kann somit in fortlaufender Simulation sehr schnell bestimmt werden. Es kann somit die Integrati-onsstufe von Partikelbeschleunigung zu Geschwindigkeit eingespart werden. Die Si-mulation der Trajektorien vereinfacht sich maßgeblich und besteht dann nur noch aus einem konvektiven Anteil der unmittelbar aus der Geschwindigkeit der Trägerströ-mung zu bestimmen ist und einer stochastisch diffusiven Störgeschwindigkeit.

4.1.3.4 Bestimmung der Zielzellen

Im Fall der Lagrangeschen Simulation bewegt sich das Partikel während eines Zeit-schritts von einem Startort im Koordinatensystem zu einem Zielort. Wesentlicher Be-standteil des Lagrangeschen Transportalgorithmus ist daher die nachträgliche Be-stimmung der Zelle, in der sich ein Partikel nach einem ausgeführten Zeitschritt befin-det. Analog zur Eulerschen Variante besteht weiterhin die Option, dass ein Partikel in der Startzelle (diejenige Zelle, in der sich das Partikel vor einem Zeitzschritt befand) verbleibt oder aber in eine andere Zelle des Rechennetzes transportiert wird.

Algorithmisch handelt es sich hierbei um eine Suchaufgabe, die über mehrere Ansät-ze erfolgen kann. Eine extensive Suche über das Netz scheidet aufgrund der hohen numerischen Komplexität der Ordnung aus.

Bei der nachbarschaftsbasierten Suche wird das Zielgebiet auf einen engen Raum um die Startzelle begrenzt. In vorliegender Entwicklung ist ein Übergang der Simulations-partikel nur bis zur ersten Nachbarschaftsebene zulässig. Abbildung 4.16 zeigt dies exemplarisch für ein 2-D unstrukturiertes Gitter. Die Zellen dieser ersten Nachbar-schaft sind dadurch gekennzeichnet, dass sie mindestens einen Gitterpunkt mit der jeweiligen Startzelle gemein haben. Aus dem begrenzten Zielgebiet leitet sich daher eine Obergrenze für den zulässigen Zeitschritt ab. Analog zu einer CFL-Bedingung beim Transport längs der Charakteristiken, ist hierbei die maximale Ausbreitungsge-schwindigkeit bedeutend:

max,

max

max

pv

rt =∆

(11)

maxt∆ ist der im Transport maximal zulässige Zeitschritt, maxr ist der kleinste Abstand

zwischen einem beliebigen Ort in der Startzelle und dem äußeren Rand der Nachbar-

zellen und max,pv ist der größtmögliche Geschwindigkeitsbetrag der Partikelbewegung.


Abbildung 4.16: Startzelle (innen) und Zielgebiet (äußere Zellen) des Lagrangeschen Transportschritts. Exemplarisch für 2-D unstrukturiertes Gitter. Grün: max. zulässiger Raumschritt. Blau: Normalenvektoren einer Zelle.

Der Test der Zellzugehörigkeit erfolgt bei vergleichbaren Algorithmen i. A. über eine eindeutige Tetraedisierung der betreffenden Zellen des Netzes (Abbildung 4.17). Wiederum der einfacheren Implementierung auf GPU-basierten Systemen geschul-det, wird eine vereinfachte Variante genutzt, die den Ort eines Partikels bezogen auf die Oberflächen der betreffenden Zelle testet. Liegt der Ort “hinter” allen Oberflächen einer Zelle, ist dies die gesuchte Zielzelle:

jNOOscal jrp ∀


Abbildung 4.17: Dekomposition einer Zelle in Tetraeder, hier exemplarisch an einer Pyramide.

Der Vereinfachung steht die erforderliche Beschränkung auf Gitter mit strikt konvexen Zellen gegenüber. Diese Einschränkung ist hinnehmbar.

4.1.3.5 Algorithmische Struktur und Auslegung der Rechenanlagen

Die Struktur der Algorithmen und die Auslegung der Rechenanlagen bedingen sich gegenseitig und können daher nicht unabhängig voneinander betrachtet werden. Es handelt sich somit um eine wechselseitige Optimierung zwischen Hardware und Soft-ware.

Aufgrund der extrem hohen Rechenleistungsanforderungen werden höchst wirtschaft-liche parallele Rechenanlangen auf Basis von GPUs (graphics processing units) ge-nutzt. Da im vorliegenden Fall zudem mehrere GPUs je Rechenknoten Einsatz finden, bestehen entsprechende Rechenanlagen aus insgesamt drei gestuften Parallelisie-rungsebenen:

• Externe Parallelisierung über Clusterknoten

• Mehrere GPUs je Rechenknoten

• Interne, massive Parallelisierung der GPUs

Abbildung 4.18 und Abbildung 4.19 zeigen den Aufbau schematisch. In dieser Anord-nung wird zum gegenwärtigen Zeitpunkt die maximale Wirtschaftlichkeit erreicht. Die Anlagen sind hierbei stark GPU-lastig. So ist die theoretische Leistungsfähigkeit der installierten GPUs ca. 100 mal höher als die der CPUs. Entsprechend muss versucht werden, die Leistung der GPUs maximal zur Entfaltung zu bringen.


Abbildung 4.18: Klassische Clusteranordnung aus Netzwerk und Rechenknoten

Abbildung 4.19: Clusterknoten bestehend je aus CPU, Host-System, mehrere (hier 3) GPUs, Netz-werkanschluss und weitere Komponenten. Typische Datentransferraten.


Jede GPU besteht intern aus einem mehrstufigen Speicher und einer Vielzahl parallel arbeitender Streamprozessoren. Die GPU stellt somit eine eigene Recheneinheit dar. Wie vorangehend bereits angedeutet, kann die extrem hohe Leistungsfähigkeit der GPUs nur sinnvoll genutzt werden, falls entsprechend einfach strukturierte, wenig granulare, rekursionsfreie Algorithmenteile vorliegen. Zudem ist eine angemessen hohe Rechenintensität erforderlich, da die interne Bandbreite der Datenanbindung zwischen Grafikspeicher und Prozessierungselementen im Vergleich zur Rechenleis-tung nur moderat ist. Grund hierfür ist der spezifisch für die ursprüngliche Verwendung als Grafikprozessoren optimierte Aufbau der GPUs. Tabelle 1 zeigt die ungefähren theoretischen Leistungen aktueller Rechnerkomponenten.

Tabelle 1: Theoretische Leistungen und Übertragungsbandbreiten der Rechnerkomponenten. Die Er-forderliche Rechenintensität gibt an, wieviele Rechenoperationen je Byte übertragener Da-ten ausgeführt werden müssten, um die Prozessierungseinheiten auslasten zu können. Ex-emplarisch für CPU: Core I-7, GPU: Nvidia GTX 580.

Speicherbandbreite (GB/s)

Rechenleistung (GF/s)

Erforderliche Intensität (F/B)

CPU (Host incl. Speicher) 12 50 ca. 4.2

GPU, residenter Speicher 190 1580 ca. 8.3

GPU über PCI 8 1580 ca. 197.5

Besonders problematisch ist die Anbindung der GPUs an das entsprechende Host-System (den betreffenden Knoten des Clusterrechners): Die Datenübertragung erfolgt über die PCI-Schnittstelle, die im Vergleich zur hohen theoretischen Leistungsfähigkeit der GPUs nur sehr niedrige Datenübertragungsraten ermöglicht. Dieser Umstand er-fordert eine Minimierung des Datenverkehrs zwischen GPUs und Host-System.

Prozessoptimierung: Aufteilung Host-System versus GPU

Zusammengefasst besteht der Ablauf der Simulationen aus folgenden, logisch abzu-grenzenden Algorithmenteilen:

1. I/O-Operationen

2. Aufbereitung der Datenstrukturen

3. Simulationsschritte zur Partikeldynamik (Koagulation und Nukleation)

4. Simulationsschritte zum Partikeltransport

5. Datenaustausch zwischen den Parallelpartitionen, wie oben.


I/O Operationen (1) müssen grundsätzlich über den betreffenden CPU-Core des Host-Systems ausgeführt werden. Die Aufbereitung der Datenstrukturen ist nur einmal zu Beginn eines Simulationslaufs erforderlich. Der Rechenaufwand ist gering, die Opera-tionen sind jedoch logisch aufwendig und hochgradig rekursiv. Eine Implementierung auf CPU-Cores ist somit selbstverständlich.

Im Gegenzug sollen numerisch intensive Algorithmenteile möglichst vollständig auf den GPU-Einheiten ausgeführt werden. Dies betrifft vor allem die stets wiederkehren-den und aufwendigen Berechnungen zur Partikeldynamik (4) in hoher Anzahl. In vor-liegender Lagrangescher Implementierung ist demgegenüber die Simulation des Par-tikeltransports (5) wesentlich weniger rechenintensiv, es besteht jedoch Datenabhän-gigkeit zu den in der Partikeldynamik erzielten Rechenfortschritten. Der Transport-schritt kann somit optional auf dem betreffenden CPU-Core des Host-Systems oder ebenfalls auf der GPU ausgeführt werden. Beide Varianten haben ausgeprägte Vor- und Nachteile. Die Auswahl erfordert zunächst eine Diskussion der Datenhaltungs-konzepte:

• A) Datenhaltung auf dem Speicher des Host-Systems

• B) Datenhaltung auf dem GPU-Speicher

• C) Hybride Datenhaltung

Zu A: Eine Datenhaltung auf dem Speicher des Host-Systems (A) böte den großen Vorteil, mit einfachsten Mitteln und geringsten Kosten ausreichende Speichermengen für eine große Partikelanzahl, etwa Millionen Zellen mit je vielen Tausend Partikeln, bereit zu stellen. In dieser Speicherstrategie würden die vergleichsweise wenig re-chenintensiven Transportschritte auf den CPU-Cores ausgeführt, da der betreffende Rechenaufwand den Datenverkehr zur GPU keinesfalls rechtfertigen würde. Leider reicht jedoch auch die Rechenintensität der Partikeldynamik bei weitem nicht aus, um die GPUs über die PCI-Schnittstelle angemessen mit Daten versorgen zu können. Eine Blockade des Programmablaufs an der PCI-Schnittstelle wäre die Folge und die Verwendung der GPUs sinnlos. Diese Variante muss daher verworfen werden.

Zu B: Werden die Daten auf den internen Grafikspeichern der GPUs selbst gehalten, kann der Datenverkehr über die PCI-Schnittstelle auf ein Minimum reduziert werden. Diesem bedeutenden Vorteil steht jedoch der Nachteil einer nur sehr begrenzten Speichermenge gegenüber: Aktuelle GPUs sind nur mit knappen Speicherausstattun-gen verfügbar, im Fall der NVidia GTX 580 sind es beispielsweise 1.5GB je Einheit. Hiervon können ca. 1.2GB für die Partikeldaten genutzt werden. Werden die Koordi-naten x,y,z und zwei weitere Eigenschaften, z. B. Partikelmasse und Sintergrad ge-speichert, sind bereits 20 Byte je Partikel erforderlich. Entsprechend können je GPU


des Typs GTX 580 die Daten für 60 Millionen Partikel aufgenommen werden. Bei 2000 Partikeln je Zelle entspricht dies einer Netzpartition mit 30.000 Zellen je GPU.

Zu C: Eine weitere, sehr interessante Option ist die hybride Datenspeicherung, d. h. die für die Partikeldynamik erforderlichen Daten werden auf der GPU abgelegt, die für einen Transportschritt erforderlichen Daten, z. B. alle Netzdaten und Partikelkoordina-ten, verbleiben auf dem Speicher des Host-Systems. In dieser Variante entfallen auf den GPUs große Anteile der je Partikel zu speichernden Datenmenge. Zudem kann ein noch etwas größerer Anteil des GPU-Speichers, im vorliegenden Fall ca. 1.4 GB statt 1.2 GB für Partikeldaten genutzt werden. Im Gegenzug wäre eine begrenzte, fort-laufende Datenübertragung über die PCI-Schnittstelle je Zeitschritt erforderlich, da einerseits die Partikeldynamik die Mobilität der entstehenden Partikel und somit die Transporteigenschaften beeinflusst, im Gegenzug der Partikeltransport über veränder-te Zellpopulationen in die Partikeldynamik zurückwirkt.

Es wird zum gegenwärtigen Zeitpunkt angenommen, dass eine hybride Datenspeiche-rung C) prinzipiell realisierbar ist. Um maximale Ausführungsgeschwindigkeiten zu ermöglichen, wurde jedoch im Rahmen des Vorhabens eine vollständig residente Speicherung auf den GPUs (B) bevorzugt. Diese Auswahl begründet auch die stark GPU-lastige Auslegung der Rechenanlage, da die Zahl der GPUs nicht nur zu einer sehr hohen Rechenleistung führt, sondern auch die maximal mögliche Größe der Si-mulationsläufe bestimmt. Im Gegenzug kann die Speicherausstattung der Hosts selbst sehr knapp gewählt werden.

Bzgl. der in vorangehendem Kapitel beschriebenen Algorithmenteile, steht eine Dis-kussion des Datenaustauschs zwischen den Parallelpartitionen (externe Parallelisie-rung) noch aus.

Externe Parallelisierung über Clusterknoten und Prozesse

Die externe Parallelisierung über die Clusterknoten erfolgt über eine klassische Ge-bietszerlegung (Englisch: domain decomposition): Das für eine Simulation erforderli-che Rechenenetz wird in Partitionen zerlegt, die je einem Prozess zugeordnet werden. Jede Partition besteht aus einem “Verantwortungsbereich” und einem darum herum-liegenden Überlappungsbereich. Letzterer überspannt typischerweise eine oder meh-rere Zellreihen, die sich ausnahmslos im Verantwortungsbereich anderer Partitionen befinden. Abbildung 4.20 zeigt eine schematische Gebietszerlegung. Die Territorien entsprechen den Verantwortungsbereichen der jeweiligen Prozesse. Abbildung 4.21 zeigt eine Detaiansicht für drei aneinander angrenzende Teilgebiete. Für das blaue Teilgebiet (links) ist ein einreihiger Überlapungsbereich skizziert. Entlang Ihrer jeweili-gen Trajektorien bewegen sich die Simulationspartikel über die Gebietsgrenzen hin-weg.


Abbildung 4.20: Skizze einer Gebietszerlegung zur parallelen Abarbeitung: Die Territorien entsprechen je dem Verantwortungsbereich eines Prozesses. Simulationspartikel entlang Ihrer Tra-jektorien über die Teilgebiete.

Abbildung 4.21: Drei aneinender angrenzende Teilgebiete (Verantwortungsgebiete: blau, rot, grün). Blaues Gebiet mit Überlappungsbereich. Simulationspartikel die in den Überlappungs-bereich eines Gebiets gelangen, werden an den jeweiligen, benachbarten Prozess weitergereicht.


Die Rechenschritte können je Partition zunächst unabhängig voneinander erfolgen, es ist jedoch in vorgegebenen Intervallen ein Datenaustausch mit benachbarten Partitio-nen erforderlich. Im Fall der Monte-Carlo Simulation führt der Partikeltransport im Randbereich zu stochastischen Übertritten in die Nachbarpartitionen. Dies ist für ein Partikel der Fall, sofern es sich nach einem Transportschritt nicht mehr im Verantwor-tungsbereich der Partition, sondern in deren Überlappungsbereich befindet.

Für jede Zelle des Überlappungsbereichs ist die betreffende Partition in dessen Ver-antwortungsbereich sich diese Zelle befindet und der dort gültige, interne Index eben-dieser Zelle bekannt. Der Datenaustausch zwischen den Partitionen wird entspre-chend nach Beendigung des Transportschritts ausgelöst, indem alle Partikel im Über-lappungsbereich an die jeweiligen Nachbarpartitionen gesendet werden.

Die Algorithmische Implementierung erfolgt unter MPI. Um eine möglichst einfache Programmstruktur zu erhalten, werden hierbei die ersten beiden Parallelisierungsebe-nen zusammengefasst, d.h. es wird jeder Parallelpartition des Netzes je ein eigener Prozess und dazugehörig je ein CPU-Core und eine GPU zugewiesen. Bzgl. der Kommunikation wird nicht explizit unterschieden, ob der für die Nachbarpartition zu-ständige Prozess auf dem gleichen oder einem anderen Clusterknoten im Netzwerk des Rechners ausgeführt wird.

Obiger Skizze, Abbildung 4.19 entsprechend, erhält jeder Clusterknoten somit je drei Partitionen des gesamten Netzes. Auf dem im Rahmen des Projekts am NST instal-lierten Hochleistungsrechner sind somit 10 Knoten je 3 GPUs verfügbar, entsprechend kann über 30 Partitionen parallelisiert werden.

Um einen Partikeltransfer zwischen zwei benachbarten Netzpartitionen ausführen zu können, müssen die betreffenden Daten auf dem Host-System des sendenden Pro-zesses verfügbar sein. Da diese entsprechend der Speicherstrategie jedoch resident auf den GPUs abgelegt sind, ist zunächst eine Übertragung von der GPU auf das Host-System erforderlich. Analog hierzu werden ankommende Partikeldaten von be-nachbarten Partitionen ebenfalls zunächst über das Host-System empfangen. Der Datenaustausch zwischen den Prozessen erfolgt also auf jeden Fall über die Host-Systeme. Dies obwohl für die Netzwerkanbindung wiederum die PCI-Schnittstelle zu-ständig ist. Der Umweg über das Host-System ist bei aktueller Technik jedoch nicht vermeidbar, bietet allerdings auch den Vorteil logisch aufwendige, jedoch wenig re-chenaufwendige, Umsortierungen auf den CPU-Cores der Host-Systeme auszuführen zu können.

Abbildung 4.22 zeigt den algorithmischen Ablauf der Datenübertragung, wie er dem schematisch skizzierten Fall, Abbildung 4.21 entspricht.


Abbildung 4.22: Datenübertragung der Simulationspartikel, entsprechend der Situation in Abbildung 4.21. Die Simulationspartikel bewegen sich ausgehend vom Gebiet des Prozesses PROC_0 zu den Nachbargebieten der Prozesse PROC_1 und PROC_2.

Basis der Programmierung

Im. Projekt sollen ausschließlich eigenständige Programme entwickelt werden, die jedoch ggf. auf offenen Paketen aufbauen. Als Ausgangspunkt stehen zur Verfügung:

▪ Das OpenFoam Paket

▪ Das Mouse Paket

▪ Das STL-Paket

▪ GPU-spezifische Programmierumgebung CUDA

▪ Offene MPI-Implementierung

Das OpenFoam Paket hat sich in den letzten Jahren stark verbreitet und wird nun-mehr für eine Vielzahl an akademischen und industriellen Problemstellungen ange-wendet, bzw. es dient als Grundlage für spezifische Weiterentwicklungen. Vorliegende Monte-Carlo Simulationsmethode nutzt die externen Datenstrukturen (File-Strukturen) der OpenFoam-Umgebung, da diese aus allen Netzgeneratoren bereitgestellt werden können. Eine Implementierung der Monte-Carlo-Methode unmittelbar innerhalb der OpenFoam-Umgebung erschien jedoch aus drei Gründen nicht sinnvoll: Erstens setzt OpenFoam in allgemeinen CFD-Anwendungen auf implizite Lösungsverfahren. Bei


diesen wird die erforderliche Rechenzeit vorwiegend für die Gleichungslösung benö-tigt. Die Diskretisierung steht somit nicht vordringlich im Fokus der Entwicklungen. Entsprechend sind die Leistungen nur moderat. Die gleiche Basis wäre im Fall der Monte-Carlo-Simulation jedoch maßgebend. Zweitens liegen kaum GPU-basierte Entwicklungen vor, so dass wesentliche Teile ohnehin neu entwickelt werden müss-ten. Drittens schließlich, ist die Programmierung von OpenFoam auf maximale Flexibi-lität ausgelegt, was sich in sehr anspruchsvollen Programmiertechniken äußert. Die OpenFoam Umgebung erscheint somit für eine extreme HPC-Entwicklung nicht ge-eignet.

Das Mouse Paket stellt ebenfalls den Anspruch, eine allgemein nutzbare Plattform für CFD-Entwicklungen zu sein. Im Gegensatz zu OpenFoam ist das Mouse Paket einfa-cher gehalten, aufgrund der geringen Verbreitung ist die Dokumentation jedoch schwach. Mouse verfügt allerdings über sehr hilfreiche Datenklassen, die das STL-Paket ergänzen und im Rahmen der vorliegenden Entwicklung ausgiebig genutzt wer-den. Eine ebenfalls unter Mouse verfügbare, sehr leistungsfähige Such- und Interpola-tionsmethode, wurde nach anfänglichen Erfolgen wieder verworfen, da die Portierbar-keit auf GPUs aufgrund der komplexen Logik nicht denkbar war.

Die Verwendung des offenen STL-Pakets (standard template library) ist selbstver-ständlich. Die GPU-spezifische Spracherweiterung CUDA ist frei verfügbar, ebenso die genutzte MPI-Implementierung.

Programmstruktur

Entsprechend den Anforderungen der Monte-Carlo Simulation auf 3-D unstrukturier-ten Gittern handelt es sich um eine extreme HPC-Anwendung, die jedoch logisch komplexe Programmteile, z. B. bezüglich der Aufbereitung der Datenstrukturen und der Datenübertragung zur Parallelisierung erfordert. Dieser vermeintliche Widerspruch wird mittels objektorientierter Programmierung unter C++ gelöst, wobei die Ausführung der numerisch intensiven Programmteile auf den GPUs lediglich von der CPU ange-stoßen werden.

Abbildung 4.23 verdeutlicht den Zusammenhang.


Abbildung 4.23: Programmiermodell. Links: Netzwerk, Mitte: logisch komplexe Programmteile die auf dem Host-System ausgeführt werden. Rechts: HPC-Teile, die durch die CPU-Cores aus der C++ Programmierung heraus angestoßen werden und auf den GPUs unter CUDA/C abgearbeitet werden.

4.1.4 Anpassung der Monte Carlo Software zur Beschreibung willkürlicher

Partikelprodukteigenschaften

Die Monte Carlo Technik wurde so implementiert, dass eine willkürliche Anzahl an Partikeleigenschaften beschrieben werden kann. Dies wird hier gezeigt anhand zweier Fallbeispiele: Koagulation und Sinterung, und Mehrkomponenten-Koagulation. Die Fallbeispiele werden hier nicht weiter untersucht oder erklärt, einziges Ziel ist, zu zei-gen, dass die Monte Carlo Software in der Lage ist, mehrere Eigenschaften gleichzei-tig zu simulieren.

Das erste Fallbeispiel ist die Beschreibung von gleichzeitiger Koagulation und Sinte-rung, anhand der Daten aus der Veröffentlichung von Kruis et al. (1993). Es wird die Änderung des Partikelvolumens aufgrund Koagulation

2

2N

dt

dN β−=

(12)

und die Änderung der Partikeloberfläche aufgrund der Sinterung simuliert:

( )

s

saa

dt

da

τ

−−=

(13)

Dieser Fall ist relevant, da damit eine 2D Verteilungsfunktion entsteht: Einmal über das Partikelvolumen, und einmal über die Partikeloberfläche (ein Indiz für den Sinter-grad). Ein Beispielergebnis ist in Abbildung 4.24 dargestellt.


Als Beleg für die Relevanz des Monte Carlo Verfahrens soll hier angeführt werden, dass das gleiche Problem in der Vergangenheit auch mittels eines sektionalen Modells mit 100 x 100 Sektionen simuliert wurde. Die Software bestand aus über 40,000 Programmzeilen, das Programm benötigte einen CRAY-Großrechner. Diese Komplexität und die Notwendigkeit, einen Großrechner zu benutzen, hatte seit der Programmerstellung verhindert, das Programm weiter einzusetzen.

Abbildung 4.24: Monte Carlo-Simulation von gleichzeitiger Koagulation und Sinterung in der Partikelentstehungsphase für einen Reaktor in dem durch SiH4 Zersetzung Si Nanopartikel entstehen. Es sind Simulationspartikel mit zwei Eigenschaften (Volumen und Sintergrad) dargestellt für verschiedene Zeiten (blau: 1 ns, rot 10 ns, grün 100 ns).

Das zweite Fallbeispiel ist die Koagulation in einem Zweikomponenten-System. Dies führt dazu, das zwei Eigenschaften notwendig sind, um die Partikel vollständig zu beschreiben: Volumen der Komponente X und Volumen der Komponente Y. Abbildung 4.25 zeigt beispielhaft die zweidimensionale Partikeleigenschaftsverteilung. Dieser Fall wurde von FST 1 bereits international veröffentlicht [68].


Abbildung 4.25: Monte Carlo-Simulation von Koagulation in einem Zwei-Komponentensystem.

4.1.5 Konzeption und Realisierung eines preiswerten Hochleistungsrechner-

systems optimiert für Monte Carlo-basierte Partikeldynamik und -

transport

Dieser Abschnitt enthält ausführliche Beschreibungen des von uns entwickelten preiswerten Hochleistungsrechnersystems, da von einigen kmUs Interesse daran be-steht. Die Anforderungen der Simulation und die Datenmengen spielen eine große Rolle bei der Auswahl der Hardware-Architektur. Die Software des Systems und das Betriebssystem spielen eine untergeordnete Rolle.

Simulationsanforderungen und Datenmengen

Im Rahmen des Projekts sollte eine Monte-Carlo Simulation die Entwicklung der rele-vanten Produkteigenschaften der Nanopartikeln beschreiben. Um die neuen Monte Carlo Simulationsmethoden einzusetzen, muss man mit großen Datenmengen und hohen Iterationszyklen rechnen.


Das System sollte in der Lage sein, Partikel-Populationen auch mit beliebigen Eigen-schaften zu berechnen. Die Form des Nanopartikel-Reaktors sollte in einem dreidi-mensionalen Simulationsgitter abgebildet werden. Es gibt keine Einschränkungen be-züglich der Gitterstruktur. Das heißt, dass Simulationsgitter könnte kartesisch, struktu-riert, unstrukturiert oder hybrid sein. Das Gitter besteht aus Elementen, die man Zel-len nennt. Die Zellen können unterschiedliche Formen haben. Diese unterschiedlichen Formen nennt man Elementtypen, wie Tetraeder, Prisma, Pyramide oder regulärer Hexaeder. Ein Millionen Zellen in dem Simulationsgitter ist das obere Ziel, sowie fünf-tausend Partikel pro Zelle.

Es müssen Annahmen getroffen werden, um die Simulationszeit innerhalb angemes-sener Grenzen zuhalten. Das Ziel ist, fünf Milliarden Partikeln in einem Reaktor zu simulieren (eine Million Zellen mit je fünftausend Partikeln). Pro Partikel braucht man 32 Byte Informationen, dreimal vier Bytes für Koordinaten, fünf mal vier Bytes für die Eigenschaften. 32 Byte für die fünfhundert Millionen Partikeln sind mehr als 14 GBy-tes. Diese Datenmenge lässt sich nur schwierig zeitlich begrenzt auf einem Rechner handhaben.

Um das numerische Konzept korrekt zu halten und die Rechenzeit zu begrenzen, sind die Partikelanzahlen auf vierhundertfünfzig Millionen reduziert worden.

Parallelisierung

Durch die physikalischen Aktionen, wie Agglomeration und Nukleation, die zwischen allen Partikel-Paaren berechnet werden müssen, steigt der Rechenaufwand schon bei wenigen Millionen simulierten Partikeln stark an. Damit übersteigen schon kleinere Szenarien selbst die Leistung moderner Rechnersysteme bei weitem. Die Rechenzeit wird nicht mehr akzeptabel. Um solche komplexen Rechnungen realisieren zu kön-nen, setzt man mehrere Prozessoren ein. In letzter Zeit geht die Entwicklung der Rechnertechnik von großen Vektorrechnern immer mehr in Richtung verteilter Spei-chersysteme (Distributed Memory Systems). Unter anderem wegen der hohen Kosten und der aufwändigen Betriebskosten (Strom und Wartung).

Die ursprüngliche Idee bestand darin, die Monte-Carlo Simulationen auf einem HPC-Workstation-Cluster laufen zu lassen. Das Simulationsgebiet sollte in etwa gleich gro-ße Gebiete zerteilt werden, die man Domains nennt. So würde man die Simulations-rechnungen parallelisiert auf mehreren Rechnern laufen lassen. Das würde die erste Parallelisierungsebene ergeben. In diesem Fall würden auf jedem Rechner 3 Domains gespeichert und 3 MPI Prozesse abgerufen. Aber, wie vorher erwähnt, wurde wegen der Datenmenge und der Iterationszahl die Lösung einer Workstation verworfen. Es wurde eine Entscheidung für einen GPU-HPC-Linux-Cluster unter Beibehaltung des


Konzeptes zur Parallelisierung getroffen. Mit der Firma Numrax wurde ein Entwurf da-für ausgelegt. Es ist ein Hybrid Linux-Cluster basierend auf einer Host-GPU Struktur und CUDA Technologie realisiert worden. Es ist ein hybrid paralleles System, bei dem die Parallelisierung auf zwei Ebenen stattfindet (Abbildung 4.26).

Abbildung 4.26: NST-Cluster Architektur.

Der Cluster wurde NSTnanoCluster genannt und besteht aus zehn Rechnern (CPU-Knoten). Jeder Rechner hat drei CUDA - fähige Grafikkarten (GPU). Nach dem Kon-zept sind die Partikeleigenschaften und die Informationen der Zellenwände (Zellenflä-chen, Normalvektor) auf den Grafikkarten zu speichern. Diese Strategie wurde aus-gewählt, um die Datentransferzeiten zwischen der CPU und GPU zu reduzieren und die Rechengeschwindigkeit zu erhöhen. Die PCIe-Bus Bandbreite zwischen CPU und GPU beträgt 8 Gigabyte pro Sekunde und die Speicherbandbreite der Grafikkarten 192.4 GB/s. Nach weiteren Abschätzungen kommt man mit 90.000 Zellen mit jeweils 5000 Partikeln zu realistischen Transferzeiten. Zum Speichern der Daten auf der Gra-fikkarte steht jeweils 1 Gigabyte freier Speicherplatz zur Verfügung.

Hardware

Die Hardware besteht aus vier wesentlichen Bestandteilen: 1. 10 Cluster-Knoten: Auf die Clusterknoten werden die parallelen Aufgabenstel-lung verteilt. Die Rechner sollen eine hohe Rechenleistung erbringen, um die Daten-mengen zu bewältigen und um die Aufgaben in kürzester Zeit zu lösen. Jeder


Clusterknoten sollte mit einem leistungsstarken Netzteil bestückt werden, um die 3 Grafikkarten ausreichend mit Strom zu versorgen. Um die Daten der Simulation zwi-schen zu speichern, ist in jedem Knoten eine Festplatte installiert (Abbildung 4.27).

Jeder Clusterknoten besteht aus: 1. Hauptplatine: Asus P6T7 WS SuperComputer Motherboard mit 7 x PCIe 2.0 x 16 slouts, mit zwei Realtek 8111C Gigabit Ethernet LAN Kontroller. 2. Arbeitsspeicher : 2 x 3072 MB DDR3. 3. CPU: Intel Core i7 quad 950 3.07 GHz. 4. Grafikkarten: 3 x GTX 580 1.5 GB. 5. Festplatte: 500 GB Hitachi SATA . 4. Netzteile: Typ Revolution 85+ 1250Watt / 1700VA ERV1250EGT.

Abbildung 4.27: Aufbau eines Cluster-Knoten. Die drei GPUs sind mit „ZOTAC“ beschriftet.

2. Server:

Der Server ist ein Rechner von dem aus die Prozesse und die benötigen Daten auf die Knoten verteilt werden. Darauf werden die benötigen Daten zentral gespeichert und visualisiert. Der Server wird auch als I/O und Verwaltungs-, Login- und Steue-rungs-Knoten eingesetzt. Ebenso als Zugangsknoten für den Fernzugriff (Remote Control) und als Firewall und Zugang zum Internet für Aktualisierung.

Er besteht aus:

1. Intel DP35DP Motherboard mit Intel 82566DC Gigabit Ethernet LAN Kontroller.

2. Intel Quad core Q6600 2.4 GHz


3. 4 x 2084MB. 4. CUDA fähige NVIDIA Geforce GT 8800 512MB Grafikkarte. 5. Zusätzlich PCI Realtek RTL8139 10/100Mbps Netzwerkkarte, mit Inter-

net verbunden. 6. DVD-R/RW Laufwerk 7. 2 x 500GB Western Digital SATA Festplatte.

3. Netzwerk:

Eine sehr zentrale Rolle im Aufbau des NSTnanoClusters spielt die Verbindungstech-nik der einzelnen Knoten untereinander und zum Server (Intrakommunikation).

Unter der Annahme der bisherigen numerischen Methoden und der Annahme, dass 5 Prozent der Partikeln aus der Zelle transportiert werden (250 Partikeln), ergeben sich (250*3000) 750.000 Partikel. In dem Worst-Case Szenario werden 50 Prozent der transportierten Partikeln die Sub-Domain verlassen. Sie werden von einem Cluster-Knoten zu einem andren Cluster-Knoten transportiert. Das sind ungefähr 375.000 Par-tikel. Damit ergibt sich ein Datenvolumen von 375.000 x 32 Byte Pro Partikel = 11,4 Megabyte.

Um die Daten zwischen den CPU-Knoten zu transportieren, wird ein Ethernet Netz-werk mit 10 Megabytes/s (100 Megabit/s) eingesetzt.

Ein Test Programm wurde geschrieben, um die MPI-Übertragungszeit in Ethernet Netzwerken zwischen zwei Cluster-Knoten zu bestimmen (Abbildung 4.28).

Abbildung 4.28: MPI-Übertragungszeit für Ethernet Netzwerke


Um die 16 Megabyte über das Ethernet Netzwerk zu transportieren, sind 1,4 Sekun-den notwendig. Im Vergleich zur Rechenzeit für die eigentliche Monte-Carlo Simulati-on ist die Netzwerktransportzeit akzeptabel. Aus diesem Grund und aus wirtschaftli-chen Gründen wurde Ethernet als Netzwerk ausgewählt.

Statt eines 100 Megabits Switches ist ein 24 Ports Gigabit Hewlett Packard (HP 24x ProCurve 2910AL-24G) Switch ausgewählt worden. So wird die Netzwerktransportzeit ungefähr auf ein Zehntel der ursprünglichen Zeit reduziert.

4. KVM Switch:

Um die zehn Cluster-Knoten bedienen und warten zu können, ist ein ATEN 16-Port KVMP Switch Model No. CS1716A eingesetzt worden. Der Server und die Cluster-Knoten können somit von einer einzigen Konsole (Tastatur, Maus und Monitor) aus gesteuert werden.

Software

Der NSTnanoCluster ist ein 64-bit Linux Cluster, auf ihm läuft ein OpenSUSE 11.3 x64 Betriebssystem. Darauf wurden mehrere Software Komponenten (tools, drivers, compilers, utilities) installiert (Abbildung 4.29). Der Cluster ist speziell entwickelt wor-den, um die nanoMC-Simulation zu realisieren und nicht als ein allgemein verwendba-res Cluster. Deshalb sind viele bekannte Cluster-Tools, wie dynamische Lastbalancie-rung oder Ressourcenverteilungen auf mehrere Aufgaben von unterschiedlichen Be-nutzern in diesem Fall nicht nötig.

Die folgen Software Komponenten wurden installiert: 1. OpenSUSE 11.3 x64 2. OpenMPI 1.5 3. GNU compiler gcc 1.4 4. SSH. 5. NVIDIA Driver 6. CUDA toolkit 4.0 7. Paraview. 8. OpenFOAM 1.7.1


Abbildung 4.29: NSTnanoCluster Software Komponenten.

Installation Das Cluster Betriebssystem ist Server Linux OpenSUSE 11.3 x64 Bit. Auf den Cluster-Knoten wurde das Betriebssystem ohne grafische Be

Autorenliste · 2012. 6. 6. · 4.3.1 Validierung der Monte Carlo Simulation der Partikeldynamik...

Documents

Transcript of Autorenliste · 2012. 6. 6. · 4.3.1 Validierung der Monte Carlo Simulation der Partikeldynamik...