Fachbereich 4: Informatik

Simulation und Visualisierung vonFluiden in einer Echtzeitanwendung mit

Hilfe der GPU

Studienarbeitim Studiengang Computervisualistik

vorgelegt von

Franz Peschel

Betreuer: Stefan Rilling(Institut für Computervisualistik, AG Computergraphik)

Koblenz, im April 2009

Erklärung

Ich versichere, dass ich die vorliegende Arbeit selbständig verfasst und keine an-deren als die angegebenen Quellen und Hilfsmittel benutzt habe.

Ja Nein

Mit der Einstellung der Arbeit in die Bibliothek bin ich einverstanden. � �

Der Veröffentlichung dieser Arbeit im Internet stimme ich zu. � �

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .(Ort, Datum) (Unterschrift)

1

Inhaltsverzeichnis

1 Einleitung 3

2 State of the Art 32.1 Das rasterbasierte Verfahren . . . . . . . . . . . . . . . . . . . . 42.2 Smoothed Particle Hydrodynamics (SPH) . . . . . . . . . . . . . 5

3 Überblick 6

4 Dynamische Fluide 64.1 Navier-Stokes-Gleichungen für inkompressible homogene Fluide . 64.2 Der Nabla-Operator . . . . . . . . . . . . . . . . . . . . . . . . . 7

4.2.1 Der Gradient . . . . . . . . . . . . . . . . . . . . . . . . 74.2.2 Die Divergenz . . . . . . . . . . . . . . . . . . . . . . . 84.2.3 Der Laplace-Operator . . . . . . . . . . . . . . . . . . . 8

4.3 Advektion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84.4 Druck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84.5 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94.6 Beschleunigung . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4.6.1 Simulation von Flüssigkeiten und Gasen . . . . . . . . . . 94.6.2 Simulation von Rauch - Auftrieb und Konvektionsströme . 94.6.3 Wirbelstärkenerhaltung . . . . . . . . . . . . . . . . . . . 10

4.7 Lösen der Navier-Stokes-Gleichungen . . . . . . . . . . . . . . . 114.7.1 Helmholtz-Hodge Zerlegung . . . . . . . . . . . . . . . . 124.7.2 Advektion . . . . . . . . . . . . . . . . . . . . . . . . . . 144.7.3 Viskose Diffusion . . . . . . . . . . . . . . . . . . . . . . 164.7.4 Beschleunigung . . . . . . . . . . . . . . . . . . . . . . . 164.7.5 Projektion - Lösen der Poisson-Gleichung . . . . . . . . . 17

4.8 Start- und Randbedingungen . . . . . . . . . . . . . . . . . . . . 18

5 Volume Rendering 195.1 Volume Slicing - texturbasiertes Volumenrendering . . . . . . . . 20

5.1.1 Object-aligned Volume Slicing . . . . . . . . . . . . . . . 205.1.2 View-aligned Volume Slicing . . . . . . . . . . . . . . . 21

6 Softwaretechnisches Konzept und Implementierung 236.1 Main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256.2 Rendermanager . . . . . . . . . . . . . . . . . . . . . . . . . . . 256.3 FluidSolverGPU3D . . . . . . . . . . . . . . . . . . . . . . . . . 25

6.3.1 init () . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256.3.2 idle_func () . . . . . . . . . . . . . . . . . . . . . . . . . 27

6.4 Implementierung . . . . . . . . . . . . . . . . . . . . . . . . . . 286.4.1 GPU - CPU . . . . . . . . . . . . . . . . . . . . . . . . . 286.4.2 Randpixel - Innenpixel . . . . . . . . . . . . . . . . . . . 29

2

6.4.3 3D Textur - Flat 3D-Textur . . . . . . . . . . . . . . . . . 296.4.4 Early-Z Culling . . . . . . . . . . . . . . . . . . . . . . . 30

6.5 Fragmentshader . . . . . . . . . . . . . . . . . . . . . . . . . . . 306.5.1 Advektion . . . . . . . . . . . . . . . . . . . . . . . . . . 306.5.2 Beschleunigung . . . . . . . . . . . . . . . . . . . . . . . 326.5.3 Temperatur . . . . . . . . . . . . . . . . . . . . . . . . . 326.5.4 Hindernisobjekte . . . . . . . . . . . . . . . . . . . . . . 346.5.5 Voxelisierung eines Objekts . . . . . . . . . . . . . . . . 346.5.6 Wirbelstärke . . . . . . . . . . . . . . . . . . . . . . . . 356.5.7 Druck und Diffusion . . . . . . . . . . . . . . . . . . . . 366.5.8 Projektion . . . . . . . . . . . . . . . . . . . . . . . . . . 37

7 Ergebnisse 38

8 Ausblick 398.1 Visualisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 398.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 408.3 Performanz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3

1 Einleitung

Mithilfe von Fluiden lassen sich einige der visuell beeindruckendsten Naturphä-nomene wie Feuer, Wasser oder Rauch realistisch darstellen. Der Simulation undVisualisierung von Fluiden wurde daher in der Vergangenheit stets großes Interes-se entgegengebracht. Fluide gelten mittlerweile in der Computergrafik als gut er-forschtes Gebiet. Fluidsimulationen wurden aber bisher eher im Offline-Renderingbeispielsweise für Animationsfilme oder auch in der Forschung und Entwicklungfür Strömungssimulationen eingesetzt. Fluidbewegungen in Echtzeit zu simulierengalt lange Zeit als zu aufwendig. Dank der jüngsten Entwicklung der Grafikhard-ware ist die Echtzeitsimulation von Fluiden nicht nur möglich, Fluidsimulationenlassen sich heute sogar in Echtzeitanwendungen wie Spiele-Engines neben anderenaufwendigen Berechnungen integrieren. Die Algorithmen wurden dabei so ange-passt, dass sie die hochparallele Architektur der GPU’s optimal ausnutzen.

Diese Arbeit liefert einen Einblick in die Technik hinter der Echtzeitsimulationvon Fluiden auf der GPU unter Verwendung von Shader Model 3.0 in OpenGL2.1. Zusätzlich wurde ein OpenGL-Renderer als Echtzeitanwendung entworfen,um die Möglichkeit der Integration der Fluidsimulation in eine solche Anwendungzu demonstrieren. Zur Ausführung der Simulation wird mindestens eine DirectX9GPU mit Shadermodel 3.0 benötigt.

2 State of the Art

Physikalische Prozesse wie die Entwicklung von Fluiden werden schon eine gan-ze Weile mithilfe von Computersystemen simuliert. Solche Simulationen habengrundsätzlich zwei verschiedene Zielstellungen. Einerseits sollen sie für wissen-schaftliche Forschung physikalische korrekt sein, andererseits sollen sie für Prä-sentationszwecke beipielsweise für den Einsatz in Computerspielen optisch plau-sibel und effizient zu berechnen sein. Eine physikalisch plausible und mit aktuellerGrafikhardware auch in Echtzeit verwendbare Simulation von Fluiden ist mit Hilfeder, für numerische Strömungssimulation geeigneten, Navier-Stokes-Gleichungenmöglich.

Die Navier-Stokes-Gleichungen wurden im 19. Jahrhundert (1827 und 1845) vonGeorge Gabriel Stokes und Claude Louis Marie Henri Navier unabhängig vonein-ander formuliert. Sie gelten auch heute noch als Grundlage der Strömungsmecha-nik und beschreiben den Impulssatz (Impulserhaltung) für newtonsche Fluide wieWasser, Gase oder Öle (viskose Flüssigkeiten) in differentieller Form [Wik09b].Zur vollständigen Beschreibung einer Fluidentwicklung verwendet man noch dieGleichung zur Massenerhaltung. Zur numerischen Lösung der Navier-Stokes-Gleichungen(ein „nichlinearen partiellen Differentialgleichungssystems 2. Ordnung“ [Wik09b])existieren unterschiedliche Ansätze.

4

Alle Elemente eines Volumens in einem Fluid mit Hilfe der Navier-Stokes-Gleichungenzu berechnen ist zu aufwendig, um praktikabel zu sein. Daher wählt man für denjeweiligen Zweck ein geeignetes Simulationsverfahren.

2.1 Das rasterbasierte Verfahren

Die Grundidee des rasterbasierten Ansatzes besteht in der Aufteilung der Szene inein gleichmäßiges 3D-Raster. Die Navier-Stokes-Gleichungen werden auf die Wer-te in den Voxelzentren des 3D-Raster angewendet. Jedes Voxelzentrum hat seineneigenen Geschwindigkeits-, Druck- und Dichtewert (Partikelanzahl). Die Dichte-werte können direkt gerendert werden. Die Werte der Volumenelemente zwischenden Voxelzentren werden interpoliert. Der Ablauf eines Simulationsschritts siehtdabei folgendermaßen aus:

• Wende die Advektion (Fluidbewegung) auf die Geschwindigkeitswerte unddie Dichtewerte an und bewege die Partikel im Fluid.

• Wende die Beschleunigung (externe Kräfte, Temperatursimulation, Wirbel-stärkenerhaltung), Hinderniserkennung und Diffusion auf die Geschwindig-keitswerte an.

• Passe den Druck im Voxel iterativ an (Inkompressibilität) [Ste08]

• Render die Dichte mittels geeignetem Volumen Rendering Verfahren.

Die rasterbasierte Methode wird hauptsächlich in der numerischen Strömungssi-mulation [Hei07] eingesetzt, da sie physikalisch plausibel ist und präzise Ergeb-nisse liefert. Die Methode ist jedoch sehr rechenaufwendig und außerdem durchdie lokale Begrenzung des 3D-Rasters nur eingeschränkt verwendbar. Der Einsatzdes Verfahrens in Computerspielen eignet sich daher am besten für Innenräume.

Abbildung 1: Hellgate London [TL07] und S.T.A.L.K.E.R.: Clear Sky [GSC08]

Dank der Entwicklung leistungsfähiger Grafikhardware im Consumerbereich wirdder rasterbasierte Ansatz bereits vereinzelt in aktuellen Computerspielen mit DirectX10-

5

Unterstützung wie z.B. dem MMORPG1 „Hellgate London“ oder dem Ego-Shooter„S.T.A.L.K.E.R.: Clear Sky“ (siehe Abb. 1) eingesetzt.

2.2 Smoothed Particle Hydrodynamics (SPH)

Die Smoothed Particle Hydrodynamics2 (kurz SPH) wurden ursprünglich von L.B. Lucy, R. A. Gingold und J. J. Monaghan für die Simulation von astrophysika-lischen Problemen entwickelt, lassen sich nach [MCG03] aber auch auf die Fluid-simulation anwenden (siehe Abb. 2). Bei den Smoothed Particle Hydrodynamicshandelt es sich um eine Lagrange-Methode, d.h. die Grundidee des Verfahrens be-steht darin, die Werte der Elemente im Volumen nicht in Voxeln zu speichern,sondern einzelne Elemente einfach im Fluid mitschwimmen zu lassen [Ste08]. EinVolumenelement in einem Fluid wird als Massepunkt aufgefasst und man wählteine Funktion, die aus dem Massepunkt wieder kleine Tropfen macht [Ste08]. Umdie Beschleunigung eines Volumenelements bestimmen zu können, werden die Po-sitionen und Beschleunigungen der Nachbarelemente verwendet. Druck entstehtdurch viele relativ nahe Nachbarelemente, Reibung entsteht durch die unterschied-lichen Geschwindigkeiten der Nachbarelemente. Zwischen den in die Berechnungeinbezogenen Elementen des Volumens kann man die Werte der übrigen Volumen-elemente einfach interpolieren [Ste08].

Abbildung 2: Smoothed Particle Hydrodynamics [MCG03]

Der partikelbasierte Ansatz ist gut für computergrafische Simulationen in Echt-zeit geeignet (z.B. Computerspiele), da man eine variable Anzahl von Partikelnbestimmen kann, die in die Berechnungen der Simulation mit einbezogen werdensollen. Bei einer geringeren Anzahl von Partikeln erreicht man eine gute Perfor-manz, allerdings auf Kosten der physikalischen Korrektheit und damit der visuel-len Qualität (O(NlogN)-Abhängigkeit des Rechenaufwands von der Teilchenzahl[Wik09c]). Das Verfahren ist außerdem laut [Wik09c] sehr dispersiv, d.h. es findet

1Massive Multiplayer Online Role Playing Game2dt.: geglättete Teilchen-Hydrodynamik [Wik09c]

6

ähnlich wie bei der Diffusion, nur stärker, ein Ausgleich zwischen den Konzentra-tionsunterschieden im Fluid statt. Dennoch kann man auch mit diesem Verfahren inKombination mit den Vortexpartikeln aus [FSWJ01] plausible visuelle Ergebnisseerzielen. Die Smoothed Particle Hydrodynamics bilden die Basis für Fluidsimula-tionen in Nvidias PhysX-Engine.

Für die Umsetzung der Aufgabenstellung der Studienarbeit fiel die Wahl auf dasrasterbasierte Verfahren trotz seiner Einschränkungen, da die Qualität der visuel-len Ergebnisse gegenüber den Ergebnissen der Smoothed Particle Hydrodynamicsüberzeugender war.

3 Überblick

Im folgenden Abschnitt 4 zu „Dynamischen Fluiden“ werden die mathematischenGrundlagen der Fluidtheorie anhand der Navier-Stokes-Gleichungen erläutert. Au-ßerdem wird die rasterbasierte numerische Lösungsmethode nach [Sta99] beschrie-ben. Im nächsten Abschnitt 5 zu „Volume Rendering“ werden dann zunächst ei-nige Techniken zur Visualisierung eines Volumendatensatzes vorgestellt. Im Ab-schnitt 6 „Softwaretechnisches Konzept und Implementierung“ wird auf Kapitel 4und 5 aufbauend die Implementierung des rasterbasierten numerischen Lösungs-verfahrens und des verwendeten Visualisierungsverfahren vorgestellt. Dabei dientein OpenGL-Renderer als Basis für die Integration der Fluidsimulation in eineEchtzeitanwendung. Die Fluidsimulation wird auf der GPU mit Hilfe von Frag-mentshadern berechnet. In den beiden letzen Abschnitten werden abschließend dieErgebnisse besprochen und ein Ausblick auf mögliche Verbesserungen und Erwei-terungen gegeben.

4 Dynamische Fluide

4.1 Navier-Stokes-Gleichungen für inkompressible homogene Fluide

Bei einem inkompressiblen homogenen Fluid3 handelt es sich um ein vereinfachtenFall des komplexen physikalischen Modells zur Simulation von Fluiden.

• Ein Fluid ist dann inkompressibel, wenn das Volumen jedes Teilbereichs desFluids konstant über die Zeit bleibt, d.h. das Volumen vergrößert sich nicht.

• Ein Fluid ist homogen, wenn seine Dichte (Partikelanzahl) d im Raum kon-stant bleibt.

• Inkompressibel und homogen bedeutet in Kombination, dass die Dichte kon-stant in Raum und Zeit bleibt.

3lat. fluidus, fließend

7

Diese Vorraussetzungen sind für dynamische Fluide üblich und schränken laut[Har03] die Verwendbarkeit des mathematischen Modells für die Simulation rea-ler Fluide wie Rauch (Gase), Wasser oder Feuer nicht ein. In den folgenden Ab-schnitten werden die Gleichungen im Detail beschrieben. Das Lösungsverfahrender Navier-Stokes-Gleichungen lässt sich in die Teilprobleme

• Advektion,

• Diffusion,

• Beschleunigung und

• Projektion

zerlegen [Har03]. Jeder dieser Teilschritte wird in den folgenden Abschnitten aus-führlich behandelt. Zunächst werden wir aber ein paar weitere Rahmenbedingungfestlegen.

Das dynamische Fluid wird in einem 3D-Raster4 mit Hilfe kartesischer Koordina-ten im Raum ~x = (x, y, z) und der Zeitvariable t simuliert. Ein Fluid mit konstanterDichte und Temperatur kann nun mit einem Geschwindigkeits-Vektor-Feld (velo-city field) ~u und einem skalaren Druckfeld (pressure field) p beschrieben werden.

~u(~x, t), p(~x, t)

Falls die Geschwindigkeit und der Druck an jedem Ort im Fluid zum Startzeitpunktt = 0 bekannt sind, kann der Zustand des Fluids zu jedem zukünftigen Zeitpunktmit Hilfe der folgenden Navier-Stokes-Gleichung für inkompressible Fluide be-stimmt werden:

∂~u

∂t= −(~u · ∇)~u− 1

d∇p + ν∇2~u+ ~F (1)

∇ · ~u = 0 (2)

Die Variable ~u und p kennen wir schon als Geschwindigkeitsfeld und Druckfeld, dist die Dichte (density), ν die Viskosität (viscosity), ~F = (fx, fy, fz) beschreibt denEinfluss externer Kräfte und ∇ = ( ∂

∂x ,∂∂y ,

∂∂z ) ist der sogenannte Nabla-Operator

(siehe Kap. 4.2). Bei der Gleichung 1 handelt es sich um die Impulserhaltungs-gleichung und bei der Gleichung 2 um die Masseerhaltungsgleichung, d.h. um dieMasseerhaltung zu gewährleisten, muss das Geschwindigkeitsfeld stets divergenz-frei (siehe Kap. 4.2) sein [Har03]. Auch die Bedingungen an den Rändern desFluids müssen berücksichtigt werden, dazu später mehr. Betrachten wir vorerst dieeinzelnen Terme der ersten Gleichung genauer.

4engl.: Grid

8

4.2 Der Nabla-Operator

In Gleichung 1 und Gleichung 2 wird der Nabla-Operator∇ auf drei verschiedeneWeisen genutzt, als Gradient-Operator, Divergenz-Operator und Laplace-Operator.

Operator Definition für den kontinuierlichen Fall

Gradient ∇p =(∂p∂x ,

∂p∂y ,

∂p∂z

)Divergenz ∇ · ~u = ∂u

∂x + ∂v∂y + ∂w

∂z

Laplace ∇2p = ∂2p∂x2 + ∂2p

∂y2 + ∂2p∂z2

Operator In Finite-Differenzen-Form für den diskreten Fall

Gradientpi+1,j,k−pi−1,j,k

2δx ,pi,j+1,k−pi,j−1,k

2δy ,pi,j,k+1−pi,j,k−1

2δz

Divergenz ui+1,j,k−ui−1,j,k

2δx ,vi,j+1,k−vi,j−1,k

2δy ,wi,j,k+1−wi,j,k−1

2δz

Laplacepi+1,j,k−2pi,j,k+pi−1,j,k

(δx)2,

pi,j+1,k−2pi,j,k+pi,j−1,k

(δy)2,

pi,j,k+1−2pi,j,k+pi,j,k−1

(δz)2

Tabelle 1: Verschieden Operatoren für die Berechnungen in der Fluidsimulation [Har05].

Die Indizes i, j, k stehen für die Koordinaten eines Voxels in einem diskretisiertendreidimensionalen Raster. δx, δy und δz bestimmen die Abstände zwischen denVoxeln im Raster in Richtung x, y und z.

4.2.1 Der Gradient

Der Gradient eines Skalarfelds ist ein Vektor, bestehend aus allen partiellen Ablei-tung des Skalarfelds [Har03].

4.2.2 Die Divergenz

Die Divergenz besitzt eine physikalische Bedeutung. Sie beschreibt, zu welchemAnteil sich Partikel in einer bestimmten Region des Raumes befinden. In den Navier-Stokes-Gleichungen wird die Divergenz auf die Geschwindigkeit des Flusses an-gewendet, um die Veränderung der Geschwindigkeit in dem Bereich um eine klei-ne Region des Fluids zu messen [Har05]. Das Fluid kann nicht zusammenge-drückt werden5. Damit diese anfangs gestellte Annahme bestand hat, wird in Glei-chung 2 sichergestellt, dass das Fluid stets Null Divergenz besitzt, d.h. Diver-genzfrei ist. Das Skalarprodukt des Divergenzoperators summiert die partiellenAbleitungen der einzelnen Komponenten des Geschwindigkeitsvektors auf. DerDivergenz-Operator kann also nur auf ein Vektorfeld wie das Geschwindigkeits-feld ~u = (u, v,w) angewendet werden [Har05].

5inkompressibel

9

4.2.3 Der Laplace-Operator

Der Gradient eines Skalarfelds ist ein Vektorfeld. Die Divergenz eines Vektorfeldsist wieder ein Skalarfeld. Wenn der Divergenz-Operator auf das Ergebnis einesGradientoperators angewendet wird, erhält man den Laplace-Operator∇·∇ = ∇2

als Ergebnis. Falls die Rastervoxel des Volumens an jeder Seite die gleiche Kan-tenlänge besitzen (kubisch), kann man den Laplace-Operator auf folgende Notationvereinfachen:

∇2p =pi+1,j,k + pi−1,j,k + pi,j+1,k + pi,j−1,k + pi,j,k+1 + pi,j,k−1 − 6pi,j,k

(δx)2

(3)

Der Laplace-Operator findet häufiger Anwendung in der Physik, beispielsweise inForm einer Diffusions-Gleichung wie der Gleichung für Hitzeentwicklung zur Be-schreibung der Hitzeverteilung in einer Region. In der Form∇2x = b wird sie auchals Poisson Gleichung bezeichnet. Im Fall b = 0 handelt es sich um die LaplaceGleichung, welche den Ursprung des Laplace-Operators darstellt. In der zweitenGleichung ∇ · ~u = 0 wird der Laplace-Operator auf ein Vektorfeld angewandt.Dabei handelt es sich um eine Vereinfachung der Notation, die Operation wirdeigentlich für jede Komponente des Vektorfelds einzeln ausgeführt [Har05].

4.3 Advektion

Der erste Term auf der rechten Seite der Gleichung 1 beschreibt die sogenannte(selbst-) Advektion6 des Geschwindigkeitsfeldes. Das Geschwindigkeitsfeld desFluids sorgt dafür, dass das Fluid Substanzen wie Partikel, Objekte oder andereDinge entlang des Strömungsflusses transportiert. Wenn man z.B. blaue Füllertintein Wasser tropft, wird die Tinte entlang des Geschwindigkeitsfeldes des Wasserstransportiert bzw. „geströmt“. Die Geschwindigkeit eines Fluids transportiert sichsogar selbst mit beim Transport der Substanz. Der beschriebene Advektions-Termsorgt also auch für die selbst-Advektion des Geschwindigkeitsfeldes [Har03].

4.4 Druck

Im zweiten Term−1d∇p wird die Beschleunigung, die durch den Druck-Gradienten

verursacht wird, beschrieben. Da Moleküle in einem Fluid sich relativ frei umein-ander herum bewegen können, prallen sie ständig zusammen oder werden einfachaneinander gedrückt. Wenn auf die Moleküle in einem lokalen Bereich eines Fluidseine Kraft ausgeübt wird, fliegen diese nicht gleich durch das gesamte Volumen desFluids, sondern prallen auf die benachbarten Moleküle, wodurch wiederum Druckensteht, da sich nun eine höhere Anzahl an Molekülen in einem Teilbereich des

6lat. advehi, heranbewegen

10

Fluids befindet. Die Moleküle werden an dieser Stelle aneinander gedrückt. Ge-nauso entsteht Unterdruck in dem Bereich, von dem die Moleküle abgezogen wer-den. Das verursacht dann ebenfalls eine Beschleunigung mit anderem Vorzeichen.Da Druck gleich Kraft pro Raumeinheit ist, führt jeglicher Druck nach Newtonszweiten Gesetz ~F = m~a zu Beschleunigung [Har03].

4.5 Diffusion

Der dritte Term ν∇2~u beschreibt die Diffusion des Fluids, deren Wirkungsgrad mitHilfe der Viskosität skaliert wird. Die Viskosität ist ein Maß dafür, wie „zähflüssig“oder dickflüssig eine Flüssigkeit fließen soll. Sirup oder Honig haben z.B. einehöhere Viskosität als Wein oder Wasser. Der so beschriebene Widerstand in einemFluid beinflusst also, wie stark sich ein Impuls in einer Flüssigkeit ausbreiten bzw.diffundieren kann.

4.6 Beschleunigung

Im vierten Term ~F wird die Beschleunigung durch externe Kräfte, die auf ein Fluidwirken, ausgedrückt. Es gibt dabei verschiedene Arten von Kräften. Globale Kräftewie die Schwerkraft unterscheiden sich von lokale Kräfte wie beispielsweise einemVentilator oder einer anderen Strömungsquelle. Damit kleine Details wie Verwir-belungen in der Strömungssimulation erhalten bleiben, wird das Fluid durch eineWirbelstärkenkraft beschleunigt (Siehe Kap. 4.6.3). Für die Simulation von Raucheignet sich zudem eine Temperatursimulation (siehe Kap. 4.6.2). Dabei wird dasFluid mit Hilfe von einer durch Wärmequellen erzeugten Auftriebskraft beschleu-nigt, die zusätzlich von der Schwerkraft, die auf die Rauchpartikel wirkt, beein-flusst wird und für sogenannte Konvektionsströmungen sorgt [Har03].

4.6.1 Simulation von Flüssigkeiten und Gasen

Die direkteste Anwendung der Fluidsimulation besteht aus der Simulation vonFlüssigkeiten und Gasen. Um eine Fluidsimulation sichtbar machen zu können,benötigt man lediglich ein Feld mit skalaren Werten, welches die zu transportie-renden Partikel enthält. Dabei kann es sich z.B. um Staub, Wassertropfen oderRußpartikel für Rauch handeln. Wir wenden den Advektions-Operator auf diesesPartikelfeld d an, um den Fluss der Partikel durch das Fluid zu simulieren.

∂~u

∂t= −(~u · ∇)~d

Dazu verwenden wir den selben Teil aus Gleichung 1, den wir für die Advektionder Geschwindigkeit nutzen. Anstelle der Geschwindigkeit ~uwird der Advektions-Operator aber auf die Dichte d angewendet [Har05].

11

4.6.2 Simulation von Rauch - Auftrieb und Konvektionsströme

Temperatur spielt eine wichtige Rolle für das Flussverhalten vieler Arten von Flui-den. Die Veränderung der Dichte und der Temperatur im Zusammenspiel erzeugenKonvektionsströme [Har03], die z.B. das Verhalten des Wetters oder auch das un-seres Kaffees in der Kaffeetasse beeinflussen. Um die dabei entstehenden Effektesimulieren zu können, müssen wir, wie in [Har03] beschrieben, den Auftrieb7 be-rechnen und ihn zu unserer Fluidsimulation hinzufügen. Dazu verwenden wir ein-fach ein neues Skalarfeld für die Temperatur T und einen Auftriebs-Operator, derdem Geschwindigkeitsfeld an der Stelle einen Impuls hinzufügt, an der die durcheine Wärmequelle erzeugte lokale Temperatur höher ist, als die Umgebungstempe-ratur Tu.

~Fbuo = σ (T− Tu)~k (4)

~k zeigt die vertikale Richtung an und σ ist ein konstanter Skalierungsfaktor. Fürdie Simulation von Rauch benötigen wir zusätzlich noch die Dichte aus dem Ska-larfeld d. Die Auftriebskraft ~Fbuo wird dabei von der Gravitationskraft, die auf dieRauchpartikel einwirkt, beeinflusst [Har03].

~Fbuo = (−µd + σ (T− Tu))~k

Bei µ handelt es sich um einen konstanten Skalierungsfaktor für die Masse derPartikel. Durch das Hinzufügen von Rauchpartikeln und Temperaturquellen (z.B.in Form einer Feuerquelle oder einer glühenden Zigarettenspitze) an eine beliebigePosition im Fluidvolumen können wir nun Rauch simulieren [Har05].

4.6.3 Wirbelstärkenerhaltung

Der Rauch einer Zigarette oder die Tinte in einem Wasserglas bilden hochdetai-lierte faszinierend anzuschauende Verwirbelungen. Durch numerischen Verlusteder Simulation verschwinden diese detailreichen Turbulenzen. Neben künstlichenAnsätzen, welche die Wirbel wieder pseudo-zufallsgesteuert hinzufügen, existierteine Ansatz von [FSWJ01] zur Wirbelstärkenerhaltung8, der die Wirbel an ihrerkorrekten Position wieder hinzufügt. Das Verfahren besitzt drei Teilschritte.

1. Im ersten Schritt wird ein Wirbelstärkefeld aus dem Geschwindigkeitsfeldberechnet. Wirbelstärke ist der Grad der Zirkulations-Rotationen in einerPartikelwolke, dargestellt mit Hilfe eines Vektorfeldes. Der Betrag einesVektors im Vektorfeld beschreibt dabei den Grad der Rotation und die Rich-tung des Vektors beschreibt die Rotationsachse der Rotation an dieser Stelle.

~Φ = ∇× ~u7engl.: buoyancy8engl.: vorticity confinement

12

2. Im zweiten Schritt werden die normalisierten Wirbelstärke-Positions-Vektorenberechnet. Sie zeigen lokal von der Region mit geringer Wirbelstärke zu derRegion mit höherer Wirbelstärke.

~X =~ω

|~ω|, ~ω = ∇|~Φ|

3. In einem dritten Schritt wird die Richtung und die Skalierung der Wirbel-stärkenkraft ~Fvc bestimmt.

~Fvc = ε(~X × ~Φ

)(5)

Wie in Abb. 3 gezeigt, kontrolliert ε den Umfang der Details, die wieder zumStrömungsfeld hinzugefügt werden [FSWJ01].

Nach der Berechnung des Wirbelstärkenkraftfelds kann es wie jede andere exter-ne Kraftquelle einfach zu dem Geschwindigkeitsfeld wieder hinzugefügt werden[Har03]. Das Hinzufügen der Wirbelstärkenerhaltung ist optional und kann ausge-lassen werden.

Abbildung 3: Hinzufügen der Wirbelstärkenerhaltung zur Simulation mit unterschiedli-chem ε-Werten. Von links nach rechts: ε = −0.25, ε = 0, ε = 0.25, ε = 0.5,ε = 0.75, ε = 1.3. Ein zu hoher Wert verursacht eine instabile Simulationund Artefakte in der Darstellung [FSWJ01].

Da sich die Qualität der Visualisieren durch das Hinzufügen der Wirbelstärkener-haltung erheblich verbessert und die Performanz nicht signifikant verschlechtert,ist das Verwenden der Wirbelstärkenerhaltung jedoch sehr zu empfehlen.

13

4.7 Lösen der Navier-Stokes-Gleichungen

Navier-Stokes-Gleichungen sind für einige vereinfachte physikalische Gegeben-heiten analytisch lösbar. Eine inkrementelle numerische Lösung erlaubt uns z.B.die Entwicklung eines Flusses über die Zeit darzustellen. Der in [Sta99] dazu vor-gestellte Algorithmus zum Lösen der Navier-Stokes-Gleichungen besteht aus meh-reren Schritten, in denen die einzelnen Gleichungen gelöst werden. Die Navier-Stokes-Gleichungen zur Impulserhaltung enthalten im dreidimensionalen Fall vierskalare Gleichungen mit vier Unbekannten: (u, v,w, p), oder einfach (~u, p) [Har05].

4.7.1 Helmholtz-Hodge Zerlegung

Ein Vektor ~v lässt sich in einzelne Basisvektoren zerlegen. Summiert man die Ba-sisvektoren auf, erhält man wieder ~v. Ein Vektor wird üblicherweise als ein n-Tupelvon Distanzen in einem n-dimensionalen Raum dargestellt. Unter Verwendung nor-malisierter Basisvektoren ~a1,~a2,~a3, ...,~an, die parallel zu den Achsen eines karte-sischen Koordinatensystems ausgerichtet sind, kann man einen Vektor als Summedarstellen: ~v = x~a1 + y~a2 + z~a3 [Har05].

Genau wie ein Vektor lässt sich auch ein Vektorfeld in eine Summe von einzelnenVektorfeldern zerlegen. Da die Geschwindigkeits- und Druckfelder zusammenhän-gen, können wir zur Vereinfachung den Druck zunächst mit Hilfe der Helmholtz-Hodge Zerlegung eliminieren, um die Anzahl der Gleichungen von vier auf drei zureduzieren [Har05]. Das Ergebniss der drei Gleichungen ist jedoch vorerst diver-gent. Laut Gleichung 2 zur Erhaltung der Masse muss das Geschwindigkeitsfeldaber divergenzfrei sein. Dieses Problem löst man wie folgt.

Der Helmholtz-Hodge Satz besagt, dass jedes Vektorfeld ~w in eine Summe vonzwei Vektorfeldern zerlegt werden kann, in ein divergenzfreies Vektorfeld ~u unddem Gradienten eines Skalarfeldes p. Die Vektoren des divergenzfreien Vektorfeldswerden dabei an den Rändern Richtung Null angenähert [Har05].

~w = ~u+∇p (6)

Das Lösen der Navier-Stokes-Gleichungen beinhaltet also nach [Sta99] für jedenZeitschritt drei Berechnungsschritte, die Advektion, Diffusion und das Hinzufü-gen eines Impulses. Das Ergebnis ist ein Geschwindigkeitsfeld ~w, dessen Diver-genz wie bereits festgestellt noch nicht Null beträgt und damit nicht die in derGleichung 2 geforderte Masseerhaltung erfüllt. Ziehen wir den Gradienten des re-sultierenden Druckfeldes in einem späteren Projektion genannten Schritt von demGeschwindigkeitsfeld ~w ab, können wir laut Helmholtz-Hodge Satz die Divergenzwieder korrigieren. Die Divergenz von ~w beträgt dann wieder Null und die Mas-seerhaltung ist gewährleistet [Har03].

14

~u = ~w −∇p

Der Helmholtz-Hodge Satz führt auch zu einer Methode, die es uns ermöglicht,das Druckfeld zu berechnen. Dazu wenden wir den Divergenz-Operator auf beidenSeiten der Gleichung ~w = ~u · p an.

∇ · ~w = ∇ · (~u+∇p)

∇ · ~w = ∇ · ~u+∇2p

Da aber die Masseerhaltungsgleichung ∇ · ~u = 0 fordert, vereinfacht sich dieGleichung zu

∇2p = ∇ · ~w (7)

mit Hilfe dieser Poisson-Gleichung lässt sich der Druck eines Fluids bestimmen.Die Gleichung wird daher auch Poisson-Druck Gleichung genannt. Nachdem wirunser divergentes Geschwindigkeitsfeld ~w bestimmt haben, können wir die Glei-chung 7 für den Druck p lösen, um mit ~w und p anschließend am Ende jedesZeitschritts wieder ein neues divergenzfreies Geschwindigkeitsfeld ~u zu berech-nen [Har05].

Als nächstes müssen wir einen Weg finden ~w zu bestimmen. Zum besseren Ver-ständnis schauen wir uns den Schritt zunächst wieder an einem Vektor an undübertragen ihn anschließend auf ein Vektorfeld. Das Skalarprodukt eines Vektor~v und eines Einheitsvektors ~e liefert uns die Projektion von ~v auf ~e. Das Skalar-produkt arbeitet also als Projektionsoperator für Vektoren. Der Operator bildet inunserem Beispiel den Vektor ~v auf seine Komponenten in Richtung des Vektors ~eab. mit Hilfe des Helmholtz-Hodge Satzes können wir einen solchen Projektions-Operator P definieren, der ein Vektorfeld ~w auf seine divergenzfreie Komponente~u projeziert. Wende wir P auf die Gleichung 6 an, erhalten wir

P ~w = P~u+ P~u+ P (∇p)

P ist laut Definition P ~w = P~u = ~u, daher ist P (∇p) = 0. Wenden wir nununseren Projektions-Operator P auf beide Seiten der Navier-Stokes-Gleichung 1zur Impulserhaltung an, erhalten wir

P∂~u

∂t= P

(−(~u · ∇)~u− 1

d∇p + ν∇2~u+ ~F

)Da ~u divergenzfrei ist, ist auch die Ableitung ∂~u

∂t divergenzfrei, wir erhalten daher

P(∂~u∂t

)= ∂~u

∂t . Der Druck-Term kann aus der Gleichung eliminiert werden, da

15

P (∇p) = 0, um die Masseerhaltung zu gewährleistet [Har05]. Wir erhalten dieGleichung

∂~u

∂t= P

(−(~u · ∇)~u+ ν∇2~u+ ~F

)(8)

Den Wert der Funktion nach der Zeit des Geschwindigkeitsfelds ~u (t) erhalten wirfür jeden Zeitpunkt t mit dem erhöhen von t um δt. Die Lösung der projezier-ten Gleichung 1 zur Impulserhaltung besteht aus mehreren Schritten [Dep08]. DasGeschwindigkeitsfeld ~u (t) zum Zeitpunkt t erhalten wir, indem wir auf der rech-ten Seite der Gleichung jeden Term einzeln auswerten und anschließend addieren.Von links nach rechts starten wir mit dem ersten Term der (selbst-) Advektion underhalten das Geschwindigkeitsfeld ~w1. ~w1 wird im zweiten Term der Diffusion ver-ändert und wir erhalten ~w2. ~w2 wird beschleunigt durch das Hinzufügen externerKräfte, Temperatur oder der Wirbelstärken-Erhaltung, wir erhalten das divergenteGeschwindigkeitsfeld ~w3 Zum Schluss wird der Projektionsoperator P auf ~w3 an-gewendet. Das resultierende Geschwindigkeitsfeld ist jetzt divergenzfrei und dieMasseerhaltung ist erfüllt. Wir lösen dazu die Gleichung 7, um den Druck p zuermitteln und den Gradienten von p von ~w3 abzuziehen [Dep08].

~u (t) Advektion−→ ~w1Diffusion−→ ~w2

Beschleunigung9−→ ~w3

Projektion−→ ~u (t + δt) [Dep08]

Jede Komponente der Gleichung 8 ist als separater Schritt zu verstehen, der einVektorfeld als Eingabe und eine verändertes Vektorfeld als Ausgabe besitzt. Diegesamte Gleichung 8 lässt sich wie P als Operator S [Har03] auffassen, der dieGleichung 8 für jeden Zeitschritt löst. Dieser Operator wird laut [Har05] als Kom-position der Operationen für

• Advektion A,

• Diffusion D,

• Beschleunigung F

• Projektion P

zu

S (~u) = P ◦ F ◦D ◦A (~u) [Har03]

Im folgenden Abschnitt werden die einzelnen Schritte einzeln erläutert und dasLösungsverfahren der Poisson-Gleichung beschrieben.

16

4.7.2 Advektion

Die Advektion beschreibt, wie ein Geschwindigkeitsfeld eines Fluids sich selbstoder andere Dinge, z.B. Partikel im Fluidvolumen transportiert. Um die Advektionvon Partikeln zu berechnen, müssen die Partikel an jeder Stelle in einem Rasteraktualisiert werden. Wir könnten die Position ~s eines jeden Partikels entlang desGeschwindigkeitsfelds über eine Zeitspanne δt bewegen.

~s (t + δt) = ~s (t) + ~uδt (9)

Dieses explizite Euler-Verfahren10 ermöglicht das explizite Integrieren von ge-wöhnlichen Differentialgleichungen und ist gleichzeitig „das einfachste und schnells-te Verfahren für das numerische Lösen eines Anfangswertproblems“11 [Wik09a].Allerdings ist es auch das ungenaueste mit der höchsten Fehlerquote. Die Diffe-rentialgleichung gehört zur Klasse der Anfangswertaufgaben, da wir sie aus vor-gegebenen Anfangsdaten ~s0 und einem Zeitpunkt t0 mit Hilfe einer Differential-gleichung ~s berechnen wollen. Natürlich gilt dabei ~s (t0) = ~s0. Neben dem Euler-Verfahren existieren auch noch genauere Verfahren wie die Runge-Kutta-Verfahrenoder die linearen Mehrschrittverfahren [Wik09a], diese sind aber etwas langsamer[Har05].

Für gewöhnliche Differentialgleichungen würden wir beispielsweise eine diskre-te Zeitspanne für die Schritte wählen. Allerdings ist der Eulersche Ansatz laut[Har05] besonders für größere Zeitschritte numerisch instabil. Der Grund für dieInstabilität des Verfahrens ist die Vorraussage der Zukunftswerte basierend auf ak-tuellen Werten, wobei in jedem Schritt Ungenauigkeitsfehler aufaddiert werden bisder Fehler zu groß wird und die Werte jeden Rahmen sprengen. Umso größer derZeitschritt gewählt wird, desto größer wird der resultierende Fehler. Fehler tretenauf, wenn der Betrag von ~u (t) δt größer ist als eine Zelle des Rasters. Ausser-dem lässt sich das Verfahren nicht auf der GPU implementieren, da die Positionender Werte in den Fragmenten nicht verändert werden können, während auf ihnengeschrieben wird. Die Partikel müssen aber für die Realisierung des Verfahren be-weglich sein [Har05].

Anstelle eines expliziten Verfahrens nutzen wir daher, wie [Sta99] vorschlägt, eininvertiertes implizites rückwärtiges Verfahren, welches uns die notwendige Stabi-lität garantiert. Im Gegensatz zum Euler-Verfahren, indem wir Partikel über einenZeitschritt advektieren, verfolgen wir jetzt die Flugbahn der Partikel jeder Raster-zelle zurück zu ihren früheren Positionen. Anschließend werden die Partikel andieser Position zu der Position ihrer Start-Rasterzelle kopiert. Dieser Rückschau-Ansatz ist laut [Sta99] für beliebige Zeitschritte und Geschwindigkeiten stabil, da

10auch Eulersches Polygonzugverfahren11auch Cauchy-Problem

17

wärend der Vorraussage der zukünftigen Werte gleichzeitig die früheren Werte kor-rigiert werden und der Gesamtfehler an einen endlichen Wert angenähert wird. Umdie Positionen der Menge q, die z.B. aus Partikeln, Geschwindigkeits- oder Tem-peraturwerten bestehen kann, im Fluid aktualisieren zu können, verwenden wirfolgende Gleichung

q (~x, t + δt) = q (~x− ~u (~x, t) δt, t) (10)

Die Gleichung beschreibt, wie wir ein Partikel entlang seines Weges zurückverfol-gen [Har05].

Abbildung 4: 2D-Advektion nach [Har03]

Abb. 4 zeigt die Advektion eines Partikels in der markierten Zelle. Das Ergebnis„x“ wird trilinear interpoliert und in die Startzelle geschrieben [Har05].

18

4.7.3 Viskose Diffusion

Viskose Fluide besitzen einen Widerstand gegen den Fluss. Dieser führt zur Diffu-sion12 der Geschwindigkeiten. Der Grad der Dämpfung der Bewegung des Flusseswird von der Viskosität ν bestimmt. Ein hoher Wert lässt das Fluid wie dickflüs-siges Öl fließen, ein niedriger Viskositätswert dagegen lässt das Fluid wie ein Gasströmen. Bei

∂~u

∂t= ν∇2~u (11)

handelt es sich um eine partielle Differentialgleichung, die die viskose Diffusionbeschreibt. Genau wie bei der Advektion verwenden wir eine Rückschau-Methode.Dazu nutzen wir die Werte des Laplace-Operator vom künftigen noch unbekanntenGeschwindigkeitsfeld. Wir gehen zunächst wieder von einer expliziten diskretenGleichung aus

~u (~x, t + δt) = ~u (~x, t) + νδt∇2~u (~xt)

Der Laplace-Operator∇2 steht hier in seiner diskreten Form. Ähnlich zur Advekti-on ist die Methode in dieser Form für große Werte von δt und ν instabil. Wir nutzendaher erneut eine in [Sta99] vorgeschlagene implizite Form der Gleichung 11:

(I − νδt∇2

)~u (~x, t+ δt) = ~u (~xt) (12)

oder auch in der Form (I − νδt∇2

)~w2 (~x) = ~w1 (~x)

I ist eine Einheitsmatrix. Die implizite Form ist wieder für beliebige Zeitschritte δtund Viskositätswerte ν stabil [Har03]. Da Diffusions-Effekte durch die Ungenau-igkeit der Berechnung im Adektionsschritt ohnehin auftreten und die Berechnungder Diffusion nur mit höheren Viskositätswerten bemerkbar ist, kann der Diffusi-onsschritt bei niedrigen Werten einfach ausgelassen werden. Die Berechnung derDiffusion ist außerdem sehr rechenaufwändig und verbessert die Qualität der Vi-sualisierung nicht, da sie Details verschwinden lässt. Wenn das Verhalten von starkviskosen Flüssigkeiten in der Simulation von Interesse ist, sollte der Diffusions-schritt verwendet werden. In allen anderen fällen ist von seiner Verwendung eherabzuraten.

4.7.4 Beschleunigung

Das Fluid bzw. das Geschwindigkeitsfeld des Fluids wird durch Hinzufügen vonexternen Kräften ~F , z.B. durch einen Ventilator, durch Auftriebskräfte ~Fbuo bei ei-ner Temperatursimulation oder durch die Wirbelstärkenkraft ~Fvc (siehe Kap. 4.6.3)

12Zerstreuung

19

beschleunigt. Das Geschwindigkeitsfeld verändert sich proportional zu den hinzu-gefügten Kräften.

~w3 = ~w2 +(~F + ~Fbuo + ~Fvc

)δt (13)

Um das Strömen eines Fluids überhaupt beobachten zu können, muss eine Kraft-quelle die Strömung erstmal erzeugen. Eine Partikelwolke sollte für die Visualisie-rung der Strömung ebenfalls hinzugefügt werden (siehe nächstes Kapitel).

4.7.5 Projektion - Lösen der Poisson-Gleichung

Der letzte Schritt der Simulation ist die Projektion P.

~u (t + δt) = P ( ~w3)

Dazu müssen zwei Poisson-Gleichungen gelöst werden:

1. Die Poisson-Druck Gleichung

2. Die Gleichung für viskose Diffusion

Wir verwenden dafür ein iteratives Verfahren, welches mit einer angenäherten Lö-sung beginnt und diese Schritt für Schritt verbessert. Die Poisson-Gleichung hatdie Form eines Matrix-Gleichungssystems

A~x = ~b

~x ist der Vektor der zu ermittelnden Werte Druck p und Geschwindigkeit ~u.~b ist einVektor von Konstanten. A ist die Matrix, die hier den Laplace-Operator ∇2 reprä-sentiert. Die iterative Lösungstechnik beginnt mit einer initialen Startschätzung x0

0

der Lösung. In jedem Schritt k wird die Lösung weiter angenähert und verbessertzu xnk . n gibt die Anzahl bereits getätigter Iterationen an. Eine einfache iterativeTechnik ist die Jacobi-Iteration. Sie ist leicht zu implementieren und daher geeig-net [Har05].

Die Gleichungen 7: ∇2p = ∇ · ~w und 12:(I − νδt∇2

)~u (~x, t+ δt) = ~u (~xt)

lassen sich beide mit Hilfe von Gleichung 3 diskretisieren und umschreiben zu derForm

~xn+1i,j,k =

~xni−1,j,k + ~xni+1,j,k + ~xni,j−1,k + ~xni,j+1,k + ~xni,j,k−1 + ~xni,j,k+1 + α~bi,j,k

β(14)

α und β sind Konstante. Die Werte für ~x,~b, α und β sind für die zwei Poisson-Gleichungen verschieden. In die Poisson-Druck Gleichung wird das Druckfeld pfür ~x und das Divergenzfeld∇ · ~w für~b eingesetzt.

20

1. Poisson-Druck Gleichung

• ~x = p

• ~b = ∇ · ~w• α = − (δx)2

• β = 6

Das Ergebnis der Poisson-Druck Gleichung ist zwar pδt statt p, der Druckwird aber nur zur Berechnung des Gradienten im Projektionsschritt genutzt.Da δt über dem gesamten Raster konstant bleibt, beeinflusst δt den Gradi-enten nicht. In die Gleichung für viskose Diffusion wird das Geschwindig-keitsfeld ~u für ~x und~b eingesetzt.

2. Gleichung für viskose Diffusion

• ~x = ~u

• ~b = ~u

• α = (δx)2

νδt

• β = 6 + α

Um die Gleichung möglichst genau zu lösen, werden eine Anzahl von Iterationendurchlaufen. Dabei wird die Gleichung 14 auf jeder Rasterzelle ausgeführt und dasErgebnis als Eingabe für den nächsten Iterationsschritt verwendet [Har03].

~x(n+1) −→ ~x(n)

Da die Jacobi-Iteration nur langsam konvergiert, müssen mindestens 20 Schritteausgeführt werden. Meine Implementation verwendet ca. 30-40 Schritte.

4.8 Start- und Randbedingungen

Für die Berechnung der Fluidentwicklung benötigen wir den Startzustand des Fluids.Als Startzustand nehmen wir das Geschwindigkeitsfeld und Druckfeld im gesam-ten Volumen als Null an. Jedes Differentialgleichungsproblem eines endlichen Be-reichs benötigt Randbedingungen. Sie werden für das Geschwindigkeitsfeld unddas Druckfeld separat definiert. Das Fluid wird in einem Raster simuliert und be-sitzt eine Grenze. Die Grenze soll in unserem Fluid durchlässig simuliert werden.Es existeren jedoch verschiedene Randbedingungen:

• no-slip BedingungWill man eine undurchlässige Begrenzung simulieren wie z.B. bei einemAquarium, nutzt man die no-slip Bedingung. Dafür wird der Wert der nächs-ten Nachbarzelle einer Randzelle mit −1 skaliert und in die Randzelle ko-piert. Die Geschwindigkeit zwischen den beiden Zellen beträgt jetzt Null.

21

• free-slip BedingungFür die Randbehandlung in einem Geschwindigkeitsfeld mit einem Hinder-nisobjekt verwendet man die free-slip Bedingung. Die Voxel in der Nähe derRänder des Hindernisobjekts nehmen die Geschwindigkeit des Hindernisob-jekts entlang seiner Oberflächennormalen an, d.h. die Oberflächennormalewird in das jeweils nächstgelegene Voxel kopiert.

• Neumann-RandbedingungFür das Druckfeld wählen wir die Neumann-Randbedingung. Der Skalie-rungsfaktor beträgt hier 1. Das führt zu einem Null-Gradienten zwischen derRandzelle und der nächsten inneren Nachbarzelle.

In dem Volumen, in dem die Rauchpartikel gespeichert werden, setzen wir in je-dem Zeitschritt am Rand den Wert Null. Dadurch können wir in beständig neueRauchpartikel hinzufügen, ohne dass das ganze Volumen nach ein paar Sekundenkomplett mit Rauchpartikeln ausgefüllt ist. Trifft ein Partikel auf den Rand, ver-schwindet es. Das Geschwindigkeitsfeld wird ebenfalls in den Randzellen stetigauf Null gesetzt. So wird der Eindruck einer unsichtbaren Barriere an den Grenzendes Fluidvolumens verhindert, an der der Rauch abzuprallen scheint.

5 Volume Rendering

Das Volumenrendern ermöglicht die Visualisierung der Strukturen im Inneren ei-nes Objekts. Es unterscheidet sich damit vom gebräuchlichen Rendern von Poly-gonoberflächen hohler Objekte ohne innerer Struktur. Ohne Berücksichtigung die-ser inneren Strukturen ist es z.B. nicht möglich, Beleuchtungsphänomene korrektdarzustellen, die nur durch diese Innenstruktur entstehen können. Besonders imBereich der medizinischen Bildverarbeitung, in dem die tatsächlichen Aufnahme-daten für die medizinische Analyse nicht manipuliert werden dürfen, benötigt manverfäschungsfreie Darstellungsverfahren. Beim Polygonrendern treten aber durchdie Triangulierung der Daten Verfälschungen auf. Daher ist dieses Verfahren nichtgeeignet. mit Hilfe des Volumenrenderverfahrens lässt sich dieses Problem lösen,denn die Verfahren verändern die Datensätze auf dem Weg zur Darstellung nicht.Neben der medizinischen Bildverarbeitung sind Volumenrenderverfahren auch füreine realistische Visualisierung von volumentrischem Nebel, Rauch oder Feuer ge-eignet.

22

Es existieren zwei verschiedene Arten von Volumenrenderverfahren

• Raycasting

• Volume Slicing

Da diese Arbeit ein Volume Slicing Verfahren verwendet, wird hier nicht näher aufdas in [CLT07] beschriebene Raycasting-Verfahren eingegangen.

5.1 Volume Slicing - texturbasiertes Volumenrendering

Beim Volume Slicing wird ein Volumen innerhalb einer Bounding Box mit Hil-fe von Slices13 dargestellt. Diese Slices sind Quad-Primitive, die zu einem Stapelgeordnet durch das Volumen in der Bounding Box gelegt werden. Die einzelnenSchichten in dem Volumen, welches beispielsweise als 3D-Textur vorliegt, werdenals 2D-Textur auf die Quads gelegt. mit Hilfe einer Alpha-Blending Funktion las-sen sich unwichtige oder durchscheinende Regionen14 ausblenden bzw. transparentmachen. Ein Volumeneffekt wird erzielt.

Beim Rendern des Stapels kommt es auf die Position der Kamera an. Bei ungünsti-gem Blickwinkel könnte man durch die einzelnen Schichten hindurchschauen. Daes beim Alpha-Blending auch auf die Reihenfolge der zu rendernden Objekte an-kommt, kann es passieren, dass die Schichten bei einer falschen Renderreihenfolgenicht mehr sichtbar sind. Ein falsche Renderreihenfolge tritt z.B. bei einer Positi-on der Kamera hinter dem Stapel auf, wenn die Schichten des Stapels von hintennach vorn gezeichnet werden. Zur Vermeidung dieser unerwünschten Effekte gibtes zwei Lösungsansätze.

5.1.1 Object-aligned Volume Slicing

Abbildung 5: Object-aligned Volume Slicing [Wil]

13Scheiben14Regionen ohne Partikel

23

Beim Object-aligned Volume Slicing (Abb. 5) wird für jede Achse ~x, ~y, ~z eineskartesischen Koordinatensystems ein Texturstapel durch das Volumen gelegt (sie-he Abb 7). Die Normalen der Ebenen eines solchen Texturstapels verlaufen par-allel zur jeweiligen Achse. Der Texturstapel, dessen Normale am ehesten parallelzur Blickrichtung der Kamera verläuft, wird gerendert. In Abb. 6 wird beispiels-weise der Texturstapel zwischen dem Schritt C und D umgeschaltet. Zusätzlichwird für korrekt funktionierendes Alpha-Blending die Reihenfolge berücksichtigt,in der die einzelnen Schichten eines Stapels gerendert werden. Befindet sich dieKamera beispielsweise hinter dem Stapel, werden die Schichten von vorn nachhinten gezeichnet. Befindet sie sich vor dem Stapel, wird dieser von hinten nachvorn gezeichnet. Das Object-aligned Volume Slicing hat den Vorteil, einfach im-plementierbar zu sein.

Abbildung 6: Je nach Kamerasicht kann zwischen den drei Texturstapeln umgeschaltetwerden [Wil]

Abbildung 7: Ein Stapel für jede Achse x,y,z

5.1.2 View-aligned Volume Slicing

Das Object-aligned Volume Slicing ist einfach zu implementieren und verhindertstörende Effekte wie z.B. das Hindurchschauen zwischen den Schichten oder dasVerschwinden des kompletten Volumes durch eine falsche Renderreihenfolge. DasHindurchschauen zwischen die Schichten wird jedoch nur verhindert, solange der

24

Öffnungswinkel15 der Kamera auf ca. 45◦-50◦ Grad eingestellt ist. Wird jedochein größerer Öffnungswinkel für eine bessere Übersicht wie z.B. 90◦ Grad ver-wendet, kann es wieder zu dem Hindurchschauen zwischen den einzelnen Schich-ten eines Stapels kommen, z.B. wenn sich das Volumen am Rand des Blickfeldesbefindet. Im Moment des Umschaltens zwischen den Stapeln entsteht beim Object-aligned Volume Slicing ausserdem ein störender visuell-wahrnehmbaren „Sprung“in der Darstellung (Poppingeffekt). Diese Effekte lassen sich mit dem in Abb. 8veranschaulichten View-aligned Volume Slicing vermeiden und die visuelle Qua-lität steigern. Allerdings ist dieses Verfahren aufwändiger zu implementieren. DieSlices werden dabei wie in Abb. 8 so durch das Volumen gelegt, dass ihre Norma-len immer parallel zur Blickrichtung der Kamera ausgerichtet sind.

Abbildung 8: View-aligned Volume Slicing [Wil]

15fov

25

6 Softwaretechnisches Konzept und Implementierung

Da die gesamte Systemarchitektur des verwendeten Renderers sehr komplex ist,werden nur die für die Fluid-Simulation wichtigen Teile des Systems beschrieben.mit Hilfe eines in Abb. 9 dargestellten Prozess-Ablaufplans der Initialisierungenund Funktionsaufrufe soll das softwaretechnische Konzept erläutert werden.

Abbildung 9: Ablaufplan eines Ausschnitts der Anwendung

26

Abbildung 10: Klassendiagramm eines Ausschnitts der Anwendung

27

6.1 Main

Startpunkt aller Prozesse ist die main.cpp. Von hier aus werden die wichtigstenObjekte und Variablen global angelegt. In der main.cpp werden auch alle notwen-digen OpenGL-Funktionen zur Initialisierung des Renderers sowie die FunktionRenderScene() aufgerufen. In RenderScene() wird in jedem Durchgangdie Szene neu berechnet, d.h. von hier aus werden die Funktionen des Renderma-nagers aufgerufen.

6.2 Rendermanager

In der Init() -Funktion des Rendermanagers, die aus der main() einmal aufge-rufen wird, werden alle Texturen, Models, die Skybox, das Terrain, einige Shaderfür die Beleuchtung, der Octree, der Fluidsolver sowie der Voxelizer geladen. DerVoxelizer erzeugt die Flat 3D-Textur für die Objekte, die sich im Fluidvolumenbefinden und von denen die Rauchpartikel abprallen sollen und übergibt sie demFluidsolver.

In jedem Durchgang der Anwendung wird aus der Funktion RenderScene()in der main.cpp die Funktion RenderGameScene() des Rendermanagers auf-gerufen. In RenderGameScene() werden die sortierten Objekte, das Terrainund das Wasser gerendert. Aus RenderScene() wird ebenfalls die FunktionRenderSmoke() des Rendermanagers aufgerufen. In RenderSmoke() wer-den die Models, die sich im Fluidvolumen befinden, gezeichnet. Ausserdem wirddie Funktion idle_func() aus dem Objekt FluidSolverGPU3D aufgerufen.idle_func() enthält dann alle weiteren notwendigen Funktionsaufrufe, um dieStrömungsentwicklung des Fluids pro Zeitschritt sowie deren Visualisierung zuberechnen.

6.3 FluidSolverGPU3D

6.3.1 init ()

In der Funktion Init() des Fluidsolvers wird nach dem Zuweisen aller Va-riablen eine Rasterstruktur GridLayer als Objekt initialisiert, die beim Setzenvon Koordinaten aus einem 3D-Datensatz in einen 2D-Datensatz benötigt wird.Der 2D-Datensatz ist eine Textur. Texturen werden ähnlich wie Arrays für CPU-Programme für GPU-Programme als Datenspeicher genutzt. Die Rasterstrukturhilft uns beim Einsortieren und Wiederfinden der 3D-Daten in einer dann als Flat3D-Textur16 (siehe Kap. 6.4.3) bezeichneten Textur. Als nächstes werden alle Tex-turen und zugehörigen Frame Buffer Objects erzeugt. Einem FBO können meh-rere Texturen zugewiesen werden. Wir benötigen hier mindestens zwei Texturenfür jedes Update eines Feldes durch ein Fragmentshader, da wir in eine Textur

16dt. flache 3D-Textur

28

nicht schreiben können, wenn sie gerade ausgelesen wird. Eine Textur wie z.B.vel_tmp dient als Zwischenspeicher. Es wäre auch möglich gewesen, für alleTexturen einen FBO zu verwenden, was Speicherplatz gespart hätte. Allerdingswären mehr Kopiervorgänge notwendig gewesen, was wiederum zu lasten der Ge-samtperformanz gegangen wäre.

Nach dem Erzeugen der Flat 3D-Textures folgt die Initialisierung der Shader, dieauf den Texturen arbeiten.

• advect3D.fragAdvektionsprogramm.

• addforce3D.fragHinzufügen externer Kräfte, z.B. ein Ventilator.

• adduniforce3D.fragHinzufügen uniformer externer Kräfte wie die Schwerkraft.

• addtemp3D.fragHinzufügen von Wärmequellen in die Temperaturtextur.

• buoyancy3D.fragBerechnen und Hinzufügen der Auftriebskraft.

• addobstacle3D.fragHinzufügen von Hindernissen im Fluidvolumen.

• vorticity3D.fragShader zur Berechnung der Wirbelstärken.

• vortforce3D.fragShader zur Berechnung der Wirbelstärkenkraft und dem Hinzufügen dieserKraft zum Geschwindigkeitsfeld.

• jacobi3D.fragBerechnet die Diffusion anhand des Viskositätsfaktors. Kann bei sehr klei-nem Viskositätsfaktor weggelassen werden. Das Programm wird auch imProjektionsschritt eingesetzt, um aus dem Divergenzfeld das Druckfeld zuberechnen.

• divergence3D.fragBerechnet das Divergenzfeld aus dem divergenten Geschwindigkeitsfeld.

• gradient3D.fragBerechnet aus dem Druckfeld und den Gradienten und zieht ihn vom diver-genten Geschwindigkeitsfeld ab. Danach ist das Geschwindigkeitsfeld diver-genzfrei.

29

• boundaryLayer3D.fragBerechnet die Werte an den Grenzen des Fluidvolumens.

• boundaryEdge3D.fragBerechnet die Werte an den äußeren Ecken des Fluidvolumens.

• volumeRender3D.fragSorgt für die Visualisierung des Fluids in einem dreidimensionalen Volumen.

30

6.3.2 idle_func ()

Die idle_func() wird in jedem Durchgang aufgerufen. Sie enthält alle Teil-schritte der Fluidberechnung. Zuerst wird der Viewport an die Auflösung der Flat3D-Textures angepasst und die Kamera so ausgerichtet, das sie orthogonal auf diegesamte Textur schaut. Jetzt können die Shader die Texturen auslesen und mit Hil-fe der FBO’s in die Texturen schreiben. Danach folgen alle Funktionen, die dieeinzelnen Berechnungsschritte der Fluidentwicklung enthalten.

Abbildung 11: Ablauf idle_func () : Funktionen und Shader

In Abb. 11 ist der Ablauf der Funktionsaufrufe in der idle_func und die in derjeweiligen Funktion aufgerufenen Shader dargestellt.

31

6.4 Implementierung

Der Renderer nutzt die Programmiersprache C++ und die Grafik-Bibliothek OpenGL2.1 mit einigen Erweiterungen17. Für die Shader wird die gut mit OpenGL kombi-nierbare und von ATI und Nvidia unterstützte Sprache GLSL verwendet.

Abbildung 12: Die einzelnen in Texturen gesepeicherten Statusfelder des Fluids. Aus-schnitte der Texturen von links nach rechts: Partikelfeld, Geschwindig-keitsfeld (grau ist keine Geschwindigkeit, grün ist Geschwindigkeit in po-sitiver y-Richtung, lila in negativer y-Richtung), Wirbelstärkenerhaltung,Druckfeld (grau ist kein Druck, weiss ist hoher Druck und schwarz niedri-ger Druck).

6.4.1 GPU - CPU

Zum Berechnen der einzelnen Entwicklungsschritte des Fluids auf der GPU wer-den wie in Abb. 12 zu sehen Texturen als Datenspeicher und Fragmentshader alsSchleifen, die auf den Pixeln der Texturen arbeiten, genutzt. Vergleichbar ist dasmit einem CPU-Programm, in denen Schleifen auf Arrayeinträgen arbeiten. Dergroße Vorteil der GPU ist dabei, dass die Pixel parallel verarbeitet werden, anstattder Reihenfolge nach, wie es bei einem CPU Programm der Fall wäre. Es könnenjedoch nur maximal soviele Pixel parallel verarbeitet werden, wie Shadereinheitendafür zu Verfügung stehen.

Ein solches Textur-Update verläuft laut [Har05] in zwei Schritten.

1. Render-to-textureEine Textur wird als Framebuffer genutzt. Die GPU kann somit direkt in dieTextur schreiben.

2. Copy-to-textureDer Framebuffer wird in eine Textur kopiert (Bei Verwendung eines FBOsfällt dieser Schritt weg)

17z.B. für FBO’s

32

Listing 1: Render-to-Texture Beispiel: Meine copy-Function12 void CFluidGPU3D::copy(TexturePtr tmp, TexturePtr x, CFBO* x_fbo, int texID)3 {4 x_fbo->Bind();5 glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT+texID);6 glClearColor(0.0, 0.0, 0.0, 0.0);7 glClear(GL_COLOR_BUFFER_BIT);89 tmp->bind(GL_TEXTURE_RECTANGLE_ARB,0);

10 glBegin(GL_QUADS);11 glMultiTexCoord2fARB(GL_TEXTURE0_ARB, 0.0f, 0.0f); glVertex2f(-1, -1);12 glMultiTexCoord2fARB(GL_TEXTURE0_ARB, fullResX, 0.0f); glVertex2f( 1, -1);13 glMultiTexCoord2fARB(GL_TEXTURE0_ARB, fullResX, fullResY); glVertex2f( 1, 1);14 glMultiTexCoord2fARB(GL_TEXTURE0_ARB, 0.0f, fullResY); glVertex2f(-1, 1);15 glEnd();1617 x_fbo->DisableFBO();18 glActiveTextureARB(GL_TEXTURE0_ARB);19 glDisable(GL_TEXTURE_RECTANGLE_ARB);20 }

Mit Hilfe der Funktion copy(...) in Listing 1 kann der Inhalt einer Textur ineine andere kopiert werden. Dabei wird Render-to-Texture mit Hilfe eines FrameBuffer Objects (kurz: FBO) eingesetzt. Ein FBO kann man anstelle des Bildschirm-Framebuffers zum Rendern in eine Textur verwenden. An einen FBO kann manmehrere Texturen, in die gerendert werden darf, binden. Mit Hilfe der Anweisun-gen in Zeile 4 bis 6 legen wir den zu verwendenden FBO und die Textur, in diegerendert werden soll, fest. Anschließend löschen wir alle Werte im Color Buf-fer des FBOs. In Zeile 9 bis 15 wird ein Viewport (Bildschirm)-füllendes Quad mitder Textur, die kopiert werden soll, gezeichnet. Man könnte jetzt noch einen Shaderauf das Quad mit der Textur anwenden. Das Quad wird direkt in die im FBO akti-vierte Textur geschrieben. In Zeile 17 wird der FBO ausgeschaltet und der Frame-buffer für den Bildschirm eingeschaltet. Nachdem der FBO ausgeschaltet wurde,kann man die Textur, in die geschrieben wurde, wieder lesen. Der Copy-to-TextureSchritt fällt weg.

6.4.2 Randpixel - Innenpixel

Da wir die Pixel am Rand (Abb. 13) einer Textur anders behandeln als die Innenpi-xel, benötigen wir für jedes Texturupdate zwei Durchgänge. Im ersten Durchgangwerden alle Innenpixel aktualisiert und im zweiten Durchgang alle Randpixel mitHilfe eines eigenen Fragmentshaders für die Randbehandlung.

6.4.3 3D Textur - Flat 3D-Textur

Ältere GPUs unterstützten zwar das Rendern von 3D-Texturen, konnten aber keineDaten in eine 3D-Textur schreiben. Daher etablierten sich die Flat 3D-Texturen.Flat 3D-Texturen sind eine praktikable und effiziente Methode zum Arbeiten mit3D-Datensätzen auf der GPU, daher werden sie in dieser Arbeit verwendet. Aller-dings unterstützen aktuelle GPUs auch das Schreiben in echte 3D-Texturen.Da wir mit einem Fluidvolumen arbeiten, müssen wir den jeweiligen Fragments-hader für jede Schicht im Volumen ausführen. Wir verwenden eine 2D-Textur,die unsere 3D-Textur repräsensiert, eine Flat 3D-Textur (Abb. 14), in der alle N

33

Abbildung 13: Randpixel - Innenpixel nach [Dep08]

Schichten des Volumens enthalten sind. Wir nutzen zum Ausführen eines Frag-mentshaders auf einer Schicht eine Schleife über alle Schichten N . N entsprichtdabei der Auflösung der Tiefe z des Volumens, d.h. bei einem würfelförmigen Vo-lumen mit Kantenlänge 64 gibt es auch 64 Schichten. Die Flat 3D-Textur wäredann (64× 8)× (64× 8) = 512× 512 Pixel.[Dep08]

Abbildung 14: 3D-Raster in einer Flat 3D-Textur nach [Dep08]

6.4.4 Early-Z Culling

mit Hilfe des Early-Z Culling Features der GPU lassen sich die Anzahl der Ope-rationen der Jacobi-Iterationen reduzieren. An Stellen im Fluid mit wenig Druckbzw. an Stellen mit relativ geringen Unterschieden der Nachbarschaftswerte kanndie Berechnung der Werte mit Hilfe der Poisson-Gleichung übersprungen werden,da sich keine Veränderung der Werte ergibt. Die Zahl der zu berechnenden Frag-

34

mente auf der GPU werden durch Ausschluss mittels des Early-Z Culling weniger.Nur die wichtigsten Regionen mit stark unterschiedlichen Nachbarschaftswertenwerden noch berechnet. Dadurch wird die Anzahl der Operationen wärend des ver-gleichsweise aufwendigen Berechnungsprozess der Poisson-Gleichung minimiert[Dep08].

6.5 Fragmentshader

6.5.1 Advektion

Das Advektionsprogramm ist der erste Fragmentshader innerhalb der idle_func(). Es transportiert die in der Textur vel gespeicherte Geschwindigkeit und spei-chert sie mit Hilfe des Framebuffer vel_fbo in Buffer-Textur vel_tmp, dieanschließend mit Hilfe der in Kapitel Render-to-texture beschriebenen Funktioncopy() wieder in die Textur vel geschrieben wird. Die Randbehandlung wirdgesondert mit Hilfe der Funktion boundary() durchgeführt. Als Eingabedatenwerden bei der Selbst-Advektion der Geschwindigkeit das Geschwindigkeitsfeldvel und als zu transportierendes Medium ebenfalls vel verwendet. Das Advekti-onsprogramm wird analog zur Selbst-Advektion noch ein zweites mal zum Trans-port der Partikel genutzt. Als Medium wird das Partikelfeld dens verwendet.

Listing 2: advect3D.frag1 uniform float dt; // Zeitschritt = 0.12 uniform float adj; // Korrekturfaktor = 0.9443 uniform float dissipation; // Verlustfaktor dissipation: 0.09974 uniform float l; // aktuelle Schicht entlang der z-Achse5 uniform vec3 resolution;6 uniform sampler2DRect velocity; // Geschwindigkeitsfeld7 uniform sampler2DRect media; // Medium8 uniform sampler2DRect offsetTexture; // Hilfstextur um die richtige Schicht9 // in z-Richtung der Flat 3D-Textur zu finden

1011 vec4 texture2DRectbilerp(in sampler2DRect media, in vec2 s)12 {13 vec4 uv;14 uv.xy = floor(s-0.5)+0.5;15 uv.zw = uv.xy + 1.0;1617 vec2 t = (s - uv.xy);1819 vec4 tex11 = texture2DRect(media, uv.xy); //(-1,-1)20 vec4 tex21 = texture2DRect(media, uv.zy); //( 1,-1)21 vec4 tex12 = texture2DRect(media, uv.xw); //(-1, 1)22 vec4 tex22 = texture2DRect(media, uv.zw); //( 1, 1)2324 // bilineare Interpolation25 return mix(mix(tex11, tex21, t.x), mix(tex12, tex22, t.x), t.y);26 }2728 vec4 advect3Drect()29 {30 //Rückverfolgen des Geschwindigkeitsfelds31 vec3 vel = -(texture2DRect(velocity, gl_TexCoord[0].xy).rgb32 - vec3(0.5, 0.5, 0.5)* adj) * dt * resolution;3334 float zValue = l + vel.z;3536 float zVal1 = floor(zValue);37 float zVal2 = zVal1 + 1.0;3839 vec2 zOffset1 = texture2DRect(offsetTexture, vec2(zVal1, 0.5)).rg;40 vec2 zOffset2 = texture2DRect(offsetTexture, vec2(zVal2, 0.5)).rg;4142 vec2 xyVal = gl_TexCoord[0].zw + vel.xy;43 xyVal = min(max(vec2(-1.0,-1.0),xyVal),resolution.xy + 1.0);44

35

45 vec4 iSample1 = texture2DRectbilerp(media, zOffset1 + xyVal);46 vec4 iSample2 = texture2DRectbilerp(media, zOffset2 + xyVal);4748 //trilineare Interpolation49 vec4 newVal = mix(iSample1, iSample2,50 (zValue/resolution.z - zVal1/resolution.z)*resolution.z);5152 return newVal;53 }5455 void main()56 {57 gl_FragColor = advect3Drect();58 }

Die Funktion advect3Drect() führt die Advektion durch. Das Geschwindig-keitsfeld und das Feld mit dem Transportmedium wird als Flat 3D-Textur velocity und media übergeben. dt gibt die Zeitspanne eines Durchgangs an. Umnumerische Genauigkeitsunterschieden zwischen verschiedenen GPU’s entgegen-zuwirken, verwende ich noch ein zur Laufzeit justierbaren zusätzlichen Parameteradj. Im Parameter l ist die aktuelle Rasterschichtnummer gespeichert. Zusammenmit der vorberechneten Textur offsetTexture, in der zu jeder Rasterschicht-nummer die zugehörigen Start-Texturkoordinaten der Flat 3D-Textur gespeichertsind, kann der Fragmentshader jeden beliebigen Voxel im Raster erreichen. DerVerlustfaktor dissipation bestimmt die Reichweite des zu transportierendenMediums bevor es verschwindet.

Abbildung 15: Trilineare Interpolation

In resolution ist die Auflösung der Flat 3D-Textur gespeichert. Da aktuelleGPU’s noch keine billineare Interpolation von floating-point Texturen unterstüt-zen, wird diese mit Hilfe der Funktion texture2DRectbilerp(...) durch-geführt. Dabei werden die vier nächstgelegenen Texel der übergebenen Textur-koordinate billinear interpoliert. Anschließend wird in advect3Drect() nochentlang der dritten Koordinate interpoliert.

6.5.2 Beschleunigung

Der Einfluss externer Kräfte findet nach der Advektion des Geschwindigkeitsfel-des statt. Ein Impuls einer externen Kraftquelle wird beispielsweise mit Hilfe des

36

Fragmentshaders in Listing 3 dem Geschwindigkeitsfeld hinzugefügt. Eine unifor-me Kraft wie die Schwerkraft kann ebenfalls Einfluss nehmen. Analog dem Hin-zufügen eines Impulses zum Geschwindigkeitsfeld werden auch neue Partikel mitHilfe des Fragmentshaders adddensity3D.frag dem Partikelfeld zugeführt.

Listing 3: addforce3D.frag12 uniform vec3 sourcepos; // Position der Kraftquelle im Fluid3 uniform vec4 sourceDir; // Staerke der Kraftquelle4 uniform vec3 windowDims;5 uniform float sourceRadius; // Radius der Kraftquelle6 uniform float l; // Aktuelle Schicht im Fluid78 uniform sampler2DRect v;9

10 float gaussian(in vec3 pos, in float radius)11 {12 return exp(-dot(pos,pos)/radius);13 }1415 vec4 splat()16 {17 vec3 pos = windowDims * sourcepos - vec3(gl_TexCoord[0].zw,l);18 float rad = windowDims.x * sourceRadius;19 vec4 c = vec4(sourceDir.rgb, clamp(dot(sourceDir,vec4(1.0,1.0,1.0,1.0)),0.0,1.0));2021 return texture2DRect(v,gl_TexCoord[0].xy) + c * gaussian(pos,rad);22 }2324 void main()25 {26 gl_FragColor = splat();27 }

6.5.3 Temperatur

Abbildung 16: Von links nach rechts: Partikelquelle, Dichte (Partikel), Wärmequelle, Ge-schwindigkeitsfeld. Im Geschwindigkeitsfeld ist der Auftrieb grün zu er-kennen. Lila eingefärbte Pixel markieren eine Abwärtsbewegung, verur-sacht durch die Gravitation, die auf die Partikel wirkt.

Wie bereits in Kap. 4.6.2 beschrieben, wird eine Temperaturtextur erstellt, in diezunächst die Wärmequellen mit Hilfe des Fragmentshaders in Listing 4 geschrie-ben werden. Anschließend wird die Auftriebskraft ~Fbuo aus Gleichung 4 im Frag-mentshader in Listing 5 mit Hilfe der Werte aus Temperaturtextur und Dichtetexturberechnet und dem Geschwindigkeitsfeld hinzugefügt (siehe Abb. 16). Die Umge-bungstemperatur und die Skalierungsfaktoren σ und µ können zur Laufzeit verän-dert werden.

37

Listing 4: addtemp3D.frag12 uniform vec3 sourcePos; // Position der Wärmequelle3 uniform float T; // Temperatur der Wärmequelle4 uniform vec3 windowDims;5 uniform float sourceRadius; // Radius der Waermequelle6 uniform float l; // Aktuelle Schicht im Fluid78 uniform sampler2DRect t; // Temperaturtextur9

10 float gaussian(in vec3 pos, in float radius)11 {12 return exp( -dot(pos, pos) / radius );13 }1415 vec4 splat()16 {17 vec3 pos = windowDims * sourcePos - vec3(gl_TexCoord[0].zw, l);18 float rad = windowDims.x * sourceRadius;19 vec4 c = vec4(T, T, T, clamp(dot(sourceTemp, vec4(1.0, 1.0, 1.0, 1.0)), 0.0, 1.0));2021 return texture2DRect(t, gl_TexCoord[0].xy) + c * gaussian(pos, rad);22 }2324 void main()25 {26 gl_FragColor = splat();27 }

Listing 5: buoyancy3D.frag12 uniform float ambientTemp; // Umgebungstemperatur3 uniform float sigma; // konstanter Skalierungsfaktor4 uniform float mu; // konstanter Skalierungsfaktor für die Partikelmasse5 uniform vec3 dir; // Zeigt die vertikale Richtung6 uniform sampler2DRect t; // Temperaturfeld7 uniform sampler2DRect d; // Dichtefeld8 uniform sampler2DRect v; // Geschwindigkeitsfeld9

10 vec4 buoyancy()11 {12 //Temperatur13 float tC = texture2DRect(t, gl_TexCoord[0].xy).r - 0.5;1415 //Dichte16 vec4 dC = texture2DRect(d, gl_TexCoord[0].xy);1718 //Geschwindigkeit19 vec4 vC = texture2DRect(v, gl_TexCoord[0].xy) - vec4(0.5,0.5,0.5,0.0);2021 //Berechnet die Auftriebskraft22 vec4 fbuo = (-mu*dC.a + sigma*(tC - ambientTemp) ) * dir;2324 //Fuegt die Auftriebskraft dem Geschwindigkeitsfeld hinzu25 vec4 u = vC + fbuo;2627 return u + vec4(0.5,0.5,0.5,0.0);28 }2930 void main()31 {32 gl_FragColor = buoyancy();33 }

6.5.4 Hindernisobjekte

Der Ablauf der Fluidsimulation verändert sich, wenn sich Hindernisobjekte imFluidvolumen befinden. Die Partikeldichte der Voxel im Hindernisobjekt beträgtNull, die Geschwindigkeitswerte dieser Voxel entsprechen der Geschwindigkeitdes Objekts. Die Bedingungen an den Rändern des Hindernisobjekts werden ge-sondert behandelt. In der Simulation wird die free-slip Bedingung verwendet. Sieverhindert, dass das Fluid in das Hindernisobjekt eindringt. Es kann stattdessen,

38

wie in Abb. 17 dargestellt, frei entlang der Oberfläche des Hindernisobjekts flie-ßen, wobei die Reibung hier ignoriert wird. Die free-slip Randbehandlung kannnach der Druckberechnung im Projektionsschritt angewandt werden. Das Ergebnisdes Projektionsschritts wird korrigiert. Wenn das aktuelle Voxel dem Rand einesHindernisobjekts angehört, verwenden wir dessen Geschwindigkeitskomponentefür das benachbarte Voxel im Fluid. Das setzt die Erkennung der Voxel, die in-nerhalb des Hindernisobjekts liegen sowie die Bestimmung ihrer Geschwindigkeitvoraus. Die Voxel innerhalb eines Hindernisobjekts werden daher vor dem Ablaufder Simulation bestimmt und in eine Flat 3D-Textur gespeichert. Diesen Prozessnennt man Voxelisierung. Die Flat 3D-Textur mit den Voxelisierungsdaten kannanschließend wie eine externe Kraft einfach dem Geschwindigkeitsfeld und demPartikelfeld hinzugefügt werden [Dep08].

Abbildung 17: Hindernisobjekt im Fluid

6.5.5 Voxelisierung eines Objekts

Für einen Teil der Voxelisierung eines Objekts (Innen-Außeninformationen) ver-wenden wir den Voxelisierungsalgorithmus aus [Dep08]. Dabei wird ein separatesRaster in Form einer RGBA Flat 3D-Textur erstellt, welches die Voxel, die sichinnerhalb eines Objekts befinden, speichert. Dafür rendert man das Objekt in einerSchleife über jede Schicht des Rasters mit Hilfe einer orthografischen Projektion.Die far clipping plane wird auf den Wert unendlich eingestellt und die near clip-

39

Abbildung 18: Die einzelnen Schichten der Voxeltextur als Ergebnis der Voxelisierungeines Autos. Die Normalen (als Farben zu sehen) werden mit Hilfe ei-nes Fragmentshaders berechnet, der einen Gradienten auf die Innen-Außeninformationen aus dem Alpha-Kanal der Voxeltextur anwendet.

ping plane auf den Tiefenwert der aktuellen Schicht. Pro Schicht wird das Objektzweimal gerendert. Zuerst wird das Objekt mit back face culling gerendert und einShader schreibt in jedes Voxel, in dem das Objekt zu sehen ist den Wert -1 in denAlpha-Kanal. Im zweiten Schritt wird das Objekt noch einmal mit front face cul-ling gerendert und mit Hilfe des Shaders wird in jedes Voxel, in dem das Objekt zusehen ist, der Wert 1 in den Alpha-Kanal geschrieben. Jetzt steht in jedem Voxel deraktuellen Schicht, das sich innerhalb des Objekts befindet, eine 1 im Alpha-Kanal.In allen übrigen steht eine 0 im Alpha-Kanal. Damit das Verfahren korrekt funk-tioniert, muss das Objekt eine geschlossene Hülle besitzen. Die Normalen und dieGeschwindigkeit eines Objekts können wir in den ersten drei Komponenten RGBder Flat 3D-Textur speichern (siehe Abb. 18). Die Geschwindigkeit eines starrenKörpers kann man über die Differenz der Vertexpositionen eines Zeitschritts ermit-teln [Dep08]. Für deformierbare Objekte benötigt man eine speziellere Behandlung(siehe dazu [CLT07]).

6.5.6 Wirbelstärke

Eine spezielle Form externer Kraftquellen ist die in [FSWJ01] beschriebene Wir-belstärkenkraft. mit Hilfe dieser Kraft können wir die kleineren detaillierten Ver-wirbelungen wieder zum Geschwindigkeitsfeld hinzufügen, die sonst durch denAdvektionsschritt verloren gehen würden. In Listing 6 ist der Fragmentshader zumBerechnen der Wirbelstärke des Geschwindigkeitsfeldes zu sehen. Listing 7 zeigtden Fragmentshader für die Berechnung der Wirbelstärkenkraft.

Listing 6: vorticity3D.frag1 uniform sampler2DRect u; //Geschwindigkeitsfeld23 vec4 vorticity3D()4 {5 vec4 L = texture2DRect(u, gl_TexCoord[0].xy - vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0);6 vec4 R = texture2DRect(u, gl_TexCoord[0].xy + vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0);7 vec4 T = texture2DRect(u, gl_TexCoord[0].xy + vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0);8 vec4 B = texture2DRect(u, gl_TexCoord[0].xy - vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0);9 vec4 U = texture2DRect(u, gl_TexCoord[1].xy) -vec4(0.5,0.5,0.5,0.0);

10 vec4 D = texture2DRect(u, gl_TexCoord[1].zw) -vec4(0.5,0.5,0.5,0.0);11 vec4 C = texture2DRect(u, gl_TexCoord[0].xy);1213 vec3 vort = vec3((T.z-B.z)-(U.y-D.y), (R.z-L.z)-(U.x-D.x), (R.y-L.y)-(T.x-B.x));1415 return vec4(vort+vec3(0.5, 0.5, 0.5), C.a);16 }

40

1718 void main()19 {20 gl_FragColor = vorticity3D();21 }

Listing 7: vortForce3D.frag1 uniform float epsilon; // Skalierung der Wirbelstärkenkraft2 uniform float timestep;3 uniform sampler2DRect vort; // Wirbelstärke4 uniform sampler2DRect u; // Geschwindigkeit56 vec4 vortForce3D()7 {8 //Auslesen der Werte Wirbelstärke in den benachbarten Voxeln9 vec4 L = abs(texture2DRect(vort, gl_TexCoord[0].xy - vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0));

10 vec4 R = abs(texture2DRect(vort, gl_TexCoord[0].xy + vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0));11 vec4 T = abs(texture2DRect(vort, gl_TexCoord[0].xy + vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0));12 vec4 B = abs(texture2DRect(vort, gl_TexCoord[0].xy - vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0));13 vec4 U = abs(texture2DRect(vort, gl_TexCoord[1].xy) -vec4(0.5,0.5,0.5,0.0));14 vec4 D = abs(texture2DRect(vort, gl_TexCoord[1].zw) -vec4(0.5,0.5,0.5,0.0));15 vec4 C = texture2DRect(vort, gl_TexCoord[0].xy) -vec4(0.5,0.5,0.5,0.0);1617 // Berechne den Gradienten des Wirbelstärkenfeldes18 vec3 gradVort = vec3(R.x-L.x, T.y-B.y, U.z-D.z);1920 vec3 V = vec3(0.0,0.0,0.0);21 vec3 vorticityForce = vec3(0.0,0.0,0.0);2223 float w = dot(gradVort,gradVort);2425 if(w!=0)26 {27 V = normalize(gradVort);28 vorticityForce = epsilon * cross(V, C.xyz);29 }3031 vec3 uNew = texture2DRect(u, gl_TexCoord[0].xy).rgb-vec3(0.5,0.5,0.5);32 uNew += timestep * vorticityForce;3334 return vec4(uNew.r+0.5, uNew.g+0.5, uNew.b+0.5, C.a);35 }3637 void main()38 {39 gl_FragColor = vortForce3D();40 }

6.5.7 Druck und Diffusion

In Kap. 4.7.5 zur Lösung der Poisson-Gleichung werden die Gleichung 7 für Druckund Gleichung 12 für Diffusion als diskretisierte Poisson-Gleichungen vorgestellt,die sich beide auch als lineares Gleichungssystem der Form A~x = ~b schreiben las-sen. Dabei ist A ein Matrix, ~x der Vektor von Unbekannten und ~b der Vektor mitden bereits bekannten Werten. Ich verwende die in Kap. 4.7.5 bereits vorgestellteeinfach Jacobi-Iteration als Technik zum Lösen der Gleichungssysteme. Ein itera-tiver Schritt zum Lösen der beiden Gleichungen für den Druck und die Diffusionkann wie in der Form der Gleichung 14

~xn+1i,j,k =

~xni−1,j,k + ~xni+1,j,k + ~xni,j−1,k + ~xni,j+1,k + ~xni,j,k−1 + ~xni,j,k+1 + α~bi,j,k

β

aus Kap. 4.7.5 beschrieben werden. In die Druckgleichung setzen wir für ~x = p · t,für~b = ∇· ~w, für α = (δx)2 und für β = 6 ein. In die Diffusionsgleichung werden

41

für ~x und~b das Geschwindigkeitsfeld ~w eingesetzt. αwird auf (δx)2

νδt und β auf 6+αgesetzt [Dep08].

Listing 8: jacobi3D.frag1 uniform float alpha;2 uniform float rBeta;3 uniform float obstacle; //Hindernisserkennung an/aus45 uniform sampler2DRect x;6 uniform sampler2DRect b;78 vec4 jacobi()9 {

10 //Auslesen der Werte der benachbarten Voxel11 vec4 xL = texture2DRect(x, gl_TexCoord[0].xy - vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0);12 vec4 xR = texture2DRect(x, gl_TexCoord[0].xy + vec2(1.0, 0.0))-vec4(0.5,0.5,0.5,0.0);13 vec4 xT = texture2DRect(x, gl_TexCoord[0].xy + vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0);14 vec4 xB = texture2DRect(x, gl_TexCoord[0].xy - vec2(0.0, 1.0))-vec4(0.5,0.5,0.5,0.0);15 vec4 xU = texture2DRect(x, gl_TexCoord[1].xy) -vec4(0.5,0.5,0.5,0.0);16 vec4 xD = texture2DRect(x, gl_TexCoord[1].zw) -vec4(0.5,0.5,0.5,0.0);1718 //Lese den Wert von b aus dem aktuellen Voxel19 vec4 bC = texture2DRect(b, gl_TexCoord[0].xy)-vec4(0.5,0.5,0.5,0.0);2021 //Schreibe Null in die Hindernisobjekt-Voxel22 if(obstacle==1){23 vec4 xC = texture2DRect(x, gl_TexCoord[0].xy)-vec4(0.5,0.5,0.5,0.0);24 if(xL.a == 1) xL = xC;25 if(xR.a == 1) xR = xC;26 if(xT.a == 1) xT = xC;27 if(xB.a == 1) xB = xC;28 if(xU.a == 1) xU = xC;29 if(xD.a == 1) xD = xC;30 }3132 //Auswerten der Jacobi-Iteration33 vec4 jac = (xL + xR + xT + xB + xU + xD + alpha * bC) * rBeta;3435 return jac + vec4(0.5, 0.5, 0.5, 0.0);36 }3738 void main()39 {40 gl_FragColor = jacobi();41 }

In Listing 8 ist der zuständige Fragmentshader für die Berechungen der Jacobi-Iteration zu sehen. Die Vektorfelder von ~x und ~b werden als Flat 3D-Texturen, αund β als uniforme floating-point Parameter übergeben. Um das skalare Druckfeldp berechnen zu können benötigen wir noch die Divergenz des Geschwindigkeits-feldes. Die Divergenz wird im Fragmentshader in Listing 9 berechnet.

Listing 9: divergence3D.frag1 uniform sampler2DRect w; // Vectorfeld, Geschwindigkeit23 vec4 divergence()4 {5 vec4 L = texture2DRect(w, gl_TexCoord[0].xy - vec2(1.0, 0.0));6 vec4 R = texture2DRect(w, gl_TexCoord[0].xy + vec2(1.0, 0.0));7 vec4 T = texture2DRect(w, gl_TexCoord[0].xy + vec2(0.0, 1.0));8 vec4 B = texture2DRect(w, gl_TexCoord[0].xy - vec2(0.0, 1.0));9 vec4 U = texture2DRect(w, gl_TexCoord[1].xy);

10 vec4 D = texture2DRect(w, gl_TexCoord[1].zw);1112 vec4 C = texture2DRect(w, gl_TexCoord[0].xy);1314 float div = (R.x - L.x + T.y - B.y + U.z - D.z) * 0.5;1516 return vec4(div+0.5, 0.5, 0.5, C.a);17 }1819 void main()20 {21 gl_FragColor = divergence();22 }

42

6.5.8 Projektion

Im Kap. 4.7.1 zur Helmholtz-Hodge Zerlegung erhält man laut Gleichung 6 ein di-vergenzfreies Geschwindigkeitsfeld, wenn man von einem divergenten Geschwin-digkeitsfeld den Gradienten eines Skalarfeldes p abzieht. Verwenden wir das Druck-feld als Skalarfeld p, können wir die geforderte Masseerhaltung gewährleisten undden Druck-Term von den Unbekannten der Navier-Stokes-Gleichung zur Impulser-haltung wie in Kap. 4.7.1 beschrieben entfernen. Listing 10 zeigt den Fragments-hader zu diesem Projektionsschritt.

Listing 10: gradient3D.frag1 uniform sampler2DRect p; // Druckfeld2 uniform sampler2DRect w; // Geschwindigkeitsfeld34 vec4 grad()5 {6 float pC = texture2DRect(p, gl_TexCoord[0].xy).r - 0.5;7 float pL = texture2DRect(p, gl_TexCoord[0].xy - vec2(1.0, 0.0)).r - 0.5;8 float pR = texture2DRect(p, gl_TexCoord[0].xy + vec2(1.0, 0.0)).r - 0.5;9 float pT = texture2DRect(p, gl_TexCoord[0].xy + vec2(0.0, 1.0)).r - 0.5;

10 float pB = texture2DRect(p, gl_TexCoord[0].xy - vec2(0.0, 1.0)).r - 0.5;11 float pU = texture2DRect(p, gl_TexCoord[1].xy).r - 0.5;12 float pD = texture2DRect(p, gl_TexCoord[1].zw).r - 0.5;1314 // Auslesen der Werte im Zentrumsvoxel und aus den benachbarten Voxeln15 // Verwendung der Werte aus den benachbarten Voxeln fuer die free-slip Randbedingung16 // an den Raendern der Hindernisobjekte17 vec4 wC = texture2DRect(w, gl_TexCoord[0].xy) - vec4(0.5,0.5,0.5,0.0);18 vec4 wL = texture2DRect(w, gl_TexCoord[0].xy - vec2(1.0, 0.0)) - vec4(0.5,0.5,0.5,0.0);19 vec4 wR = texture2DRect(w, gl_TexCoord[0].xy + vec2(1.0, 0.0)) - vec4(0.5,0.5,0.5,0.0);20 vec4 wT = texture2DRect(w, gl_TexCoord[0].xy + vec2(0.0, 1.0)) - vec4(0.5,0.5,0.5,0.0);21 vec4 wB = texture2DRect(w, gl_TexCoord[0].xy - vec2(0.0, 1.0)) - vec4(0.5,0.5,0.5,0.0);22 vec4 wU = texture2DRect(w, gl_TexCoord[1].xy) - vec4(0.5,0.5,0.5,0.0);23 vec4 wD = texture2DRect(w, gl_TexCoord[1].zw) - vec4(0.5,0.5,0.5,0.0);2425 // Geschwindigkeitsfeld und Maske des Geschwindigkeitsfelds durch26 // die free-slip Randbedingung27 vec3 obstV = vec3(0,0,0);28 vec3 vMask = vec3(1,1,1);2930 // 1. Elimination des unkorrekten Druck-Gradient Terms31 // (d.h. Druck Werte im Hindernisobjekt werden ignoriert)32 // 2. berechne Geschwindigkeit durch free-slip Randbedingung33 if(wL.a == 1){pL = pC; obstV.x = wL.x; vMask.x = 0;}34 if(wR.a == 1){pR = pC; obstV.x = wR.x; vMask.x = 0;}35 if(wB.a == 1){pB = pC; obstV.y = wB.y; vMask.y = 0;}36 if(wT.a == 1){pT = pC; obstV.y = wT.y; vMask.y = 0;}37 if(wU.a == 1){pU = pC; obstV.z = wU.z; vMask.z = 0;}38 if(wD.a == 1){pD = pC; obstV.z = wD.z; vMask.z = 0;}3940 // Gradient41 vec3 gradP = vec3(pR - pL, pT - pB, pU - pD);4243 // Projektion44 vec4 u = vec4(0.0,0.0,0.0,0.0);45 u.xyz = wC.xyz - gradP;4647 // Geschwindigkeitskorrektur durch free-slip Randbedingung48 u.xyz = (vMask * u.xyz) + obstV;4950 return u + vec4(0.5,0.5,0.5,0.0);51 }5253 void main()54 {55 gl_FragColor = grad();56 }

In das Fluidvolumen können Hindernisse in Form von Objekten platziert werden.Die Bedingungen an den Rändern eines solchen Objekts werden gesondert behan-delt. Für das Druckfeld verwenden wir die Neumann-Randbedingung. Der Gradi-ent des Druckfeldes entlang der Oberflächennormalen des Hindernisobjekts sollte

43

dabei den Wert Null annehmen. Beim Geschwindigkeitsfeld verwenden wir ei-ne free-slip Randbedingung, damit die Geschwindigkeit des Objekts auf das Ge-schwindigkeitsfeld übertragen werden kann [Dep08].

7 Ergebnisse

In der Studienarbeit konnte gezeigt werden, dass die physikalisch basierte Simu-lation von dreidimensionalen Fluiden in Echtzeit auf aktuellen GPU’s möglich ist.Die Fluidsimulation wurde erfolgreich in eine Echtzeitanwendung integriert undkann mit Objekten interagieren. Die Echtzeitanwendung ist ein kleiner OpenGL-Renderer, mit dessen Hilfe man Terrain, Objekte, Wasser usw. rendern kann. DieAnwendung erlaubt die Manipulation der Parameter der Fluidsimulation und dasSpeichern der aktuellen Konfiguration. Der Nutzer kann sich mit den vielfältigenMöglichkeiten und faszinierenden Effekten einer Fluidsimulation vertraut machen.Je nach GPU kann die Rasterauflösung des Fluids angepasst werden. Auf dem Test-system mit einem Grafikbeschleuniger von Nvidia (GeForce GTX260) erreichtedie Simulation bei einer Rasterauflösung von 100×100×100 Pixeln durchschnitt-lich 31 fps.

Die folgenden Abbildungen sind Screenshots aus der eigens erstellten Echtzeitan-wendung in Form eines OpenGL-Renderers, in die die Fluidsimulation integriertwurde. Sie demonstrieren einige Effekte wie Explosionen (Abb. 20), Flammen(Abb. 21), Nebel (Abb. 25, Abb. 24, Abb. 23), Rauchschwaden (Abb. 22, Abb.19, Abb. 26, Abb. 27) oder das Verhalten des Fluids an Hindernisobjekten (Abb.28, Abb. 29) die mit einer Fluidsimulation visualisiert werden können.

8 Ausblick

Die vorliegende Fluidsimulation bietet vielfältige Erweiterungsmöglichkeiten.

8.1 Visualisierung

Für eine verbesserte Volumenvisualisierung könnte beispielsweise das Raycasting-Verfahren aus [CLT07] implementiert werden. Eine größere Herausforderung wäreauch die Suche nach einem geeigneten Beleuchtungsverfahren. Das Kapitel „CloudIllumination and Rendering“ aus der Dissertation von [Har03] liefert dazu ausführ-lichere Informationen. Beim View-aligned Volume Slicing sind an den Stellen, andenen die Slices die Geometrie der Szene schneiden, Artefakte, wie z.B. unnatür-lich scharfe Kanten, zu sehen. Zur Vermeidung dieser Artefakte könnte man die„Soft Particles“-Methode aus [Lor] verwenden.

44

Abbildung 19: Visualisierung eines Vulkanausbruchs

8.2 Simulation

Neben der Simulation von Rauch- und Gaseffekten kann die Fluidsimulation auchals Basis für die Realisierung von Feuer oder Wassereffekten verwendet werden[CLT07]. Im Beschleunigungsschritt kann man zudem eine Temperatursimulation[Har03] hinzufügen. Nach der Beschleunigung und vor der Projektion kann für dieSimulation von Dynamischen Wolken ein Thermodynamik-Schritt wie in [Har03]beschrieben eingefügt werden. Für deformierbare Objekte ließe sich ein verbes-sertes Voxelisierungsverfahren [CLT07] unter Verwendung des Geometry Shadersrealisieren. Der Renderer der Applikation bietet zusätzlich zahlreiche Experimen-tiermöglichkeiten für den Einsatz von Fluiden in Echtzeitanwendungen.

8.3 Performanz

Im Projektionsschritt lassen sich durch eine Vektorisierung der skalaren Druckwer-te kleinere Texturen verwenden. Dazu werden immer vier skalare Werte in einenRGBA-Vektor gepackt. So erhält man weniger Fragmente, die der Shader wärendder Jacobi Iteration berechnen muss. Mehr dazu in [Har03].

45

Abbildung 20: Visualisierung einer Explosion

Abbildung 21: Visualisierung einer bunten Flamme

Abbildung 22: Visualisierung von Rauchschwaden

Abbildung 23: Verändert man den Blendmode, erscheinen die Rauchschwaden wieDampf

Abbildung 24: Visualisierung bläulich schimmernder Nebelschwaden

Abbildung 25: Farbige Nebelschwaden

Abbildung 26: Die Entwicklung eines Rauchpilzes in einzelnen Schritten in einem 150×150× 150 großen Grid.

Abbildung 27: Visualisierung eines Rauchpilzes

Abbildung 28: Hinderniserkennung in einzelnen Zeitschritten

Abbildung 29: Hindernisobjekt im Fluid

Abbildungsverzeichnis

1 Hellgate London [TL07] und S.T.A.L.K.E.R.: Clear Sky [GSC08] 42 Smoothed Particle Hydrodynamics [MCG03] . . . . . . . . . . . 53 Wirbelstärkenerhaltung mit unterschiedlichen ε-Werten . . . . . . 114 2D-Advektion nach [Har03] . . . . . . . . . . . . . . . . . . . . 155 Object-aligned Volume Slicing [Wil] . . . . . . . . . . . . . . . . 206 Texturstapel umschalten . . . . . . . . . . . . . . . . . . . . . . 217 Ein Stapel für jede Achse x,y,z . . . . . . . . . . . . . . . . . . . 218 View-aligned Volume Slicing [Wil] . . . . . . . . . . . . . . . . . 229 Ablaufplan eines Ausschnitts der Anwendung . . . . . . . . . . . 2310 Klassendiagramm eines Ausschnitts der Anwendung . . . . . . . 2411 Ablauf idle_func () : Funktionen und Shader . . . . . . . . . . . . 2712 Statusfelder der Fluidsimulation . . . . . . . . . . . . . . . . . . 2813 Eandpixel - Innenpixel nach [Dep08] . . . . . . . . . . . . . . . . 2914 3D-Raster in einer Flat 3D-Textur nach [Dep08] . . . . . . . . . . 3015 Trilineare Interpolation . . . . . . . . . . . . . . . . . . . . . . . 3216 Temperatur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3317 Hindernisobjekt im Fluid . . . . . . . . . . . . . . . . . . . . . . 3418 Voxeltextur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3519 Visualisierung eines Vulkanausbruchs . . . . . . . . . . . . . . . 3920 Visualisierung einer Explosion . . . . . . . . . . . . . . . . . . . 4121 Visualisierung einer bunten Flamme . . . . . . . . . . . . . . . . 4122 Visualisierung von Rauchschwaden . . . . . . . . . . . . . . . . 4223 Dampf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4224 Visualisierung bläulich schimmernder Nebelschwaden . . . . . . 4325 Farbige Nebelschwaden . . . . . . . . . . . . . . . . . . . . . . . 4326 Entwicklung eines Rauchpilzes . . . . . . . . . . . . . . . . . . . 4427 Visualisierung eines Rauchpilzes . . . . . . . . . . . . . . . . . . 4428 Hinderniserkennung in einzelnen Zeitschritten . . . . . . . . . . . 4529 Hindernisobjekt im Fluid . . . . . . . . . . . . . . . . . . . . . . 45

51

Listings

1 Render-to-Texture Beispiel: Meine copy-Function . . . . . . . . . 292 advect3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303 addforce3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . 324 addtemp3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . 335 buoyancy3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . 336 vorticity3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . 357 vortForce3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . 368 jacobi3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . . 369 divergence3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . 3710 gradient3D.frag . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

52

Literatur

[CLT07] CRANE, Keenan ; LLAMAS, Ignacio ; TARIQ, Sarah: Real-Time Simu-lation and Rendering of 3DFluids. In: NGUYEN, Hubert (Hrsg.): GPUGems 3. Addison Wesley Professional, August 2007, Kapitel 30

[Dep08] DEPARTMENT OF CONTROL ENGINEERING AND INFORMATION

TECHNOLOGY: Gridfluid - 3D Grid Based Fluid Simulator for Hewlett-Packard Scalable Visualization Array / Budapest University of Techno-logy and Economics. 2008. – Forschungsbericht

[FSWJ01] FEDKIW, Ronald ; STAM, Jos ; WANN JENSEN, Henrik: Visual Si-mulation of Smoke. In: Proceedings of ACM SIGGRAPH 2001, 2001(Computer Graphics Proceedings, Annual Conference Series), S. 15–22

[GSC08] GSC GAME WORLD: S.T.A.L.K.E.R.: Clear Sky. http://stalker.deepsilver.com. Version: 2008. – S.T.A.L.K.E.R.: Clear Sky HDDirectX10 Feature Demo Video

[Har03] HARRIS, Mark J.: Real-Time cloud simulation and rendering, Universityof North Carolina at Chapel Hill, Diss., 2003

[Har05] HARRIS, Mark J.: Fast fluid dynamics simulation on the GPU. In: SIG-GRAPH ’05: ACM SIGGRAPH 2005 Courses, ACM Press, 2005

[Hei07] HEINECKE, Berthold: Physikalische Rauchsimulation auf Partikelbasisin Echtzeit mit der PhysX-Engine. September 2007. – Belegarbeit an derTechnischen Universität Dresden

[Lor] LORACH, Tristan: Soft Particles. – NVIDIA DirectX 10 SDK

[MCG03] MÜLLER, Matthias ; CHARYPAR, David ; GROSS, Markus: Particle-based fluid simulation for interactive applications. In: SCA ’03: Procee-dings of the 2003 ACM SIGGRAPH/Eurographics symposium on Com-puter animation. Aire-la-Ville, Switzerland, Switzerland : EurographicsAssociation, 2003. – ISBN 1–58113–659–5, S. 154–159

[Sta99] STAM, Jos: Stable Fluids. In: ROCKWOOD, Alyn (Hrsg.): Siggraph1999, Computer Graphics Proceedings. Los Angeles : Addison WesleyLongman, 1999, S. 121–128

[Ste08] STELZMANN, Robert: PPartikelsysteme und Fluidsimulationen. Au-gust 2008. – Proseminar Computergrafik an der Technischen UniversitätDresden

[TL07] TARIQ, Sarah ; LLAMAS, Ignacio: Real-Time Volumetric Smoke usingD3D10. Proceedings of the Game Developer Conference 2007, März2007. – NVIDIA Developer Technology

53

[Wik09a] WIKIPEDIA: Eulersches Polygonzugverfahren— Wikipedia, Die freieEnzyklopädie. http://de.wikipedia.org/w/index.php?tit\le=Eulersches_Polygonzugverfahren\&oldid=58725467. Version: 2009. – [Online; Stand 6. April 2009]

[Wik09b] WIKIPEDIA: Navier-Stokes-Gleichungen— Wikipedia, Die freieEnzyklopädie. http://de.wikipedia.org/w/index.php?title=\Navier-Stokes-Gleichungen\&oldid=58713318. Version: 2009. – [Online; Stand 5. April 2009]

[Wik09c] WIKIPEDIA: Smoothed Particle Hydrodynamics— Wikipedia, Die freieEnzyklopädie. http://de.wikipedia.org/w/index.php?title=Smooth\ed_Particle_Hydrodynamics\&oldid=56812835. Version: 2009. – [Online; Stand 6. April 2009]

[Wil] WILLKOMM, Dennis: Visuelle Effekte mit volumetrischen Shadern. –Studienarbeit an der Universität Koblenz-Landau

54