Partikelbasierte Oberflächenerosionssimulation


Hinweis: Der vollständige Quellcode des Projekts sowie Erläuterungen zu seiner Verwendung und zum Lesen finden Sie auf Github [ hier ].

Ich machte eine Pause von meiner Masterarbeit, um an dem zu arbeiten, was ich lange aufgeschoben hatte: Verbesserte Geländegenerierung für mein Territory-Projekt . Eine einfache Möglichkeit, es zu implementieren, ist die hydraulische Erosion, weshalb ich es erstellt habe!

Für ein eintägiges Software-Puzzle funktionierte es ziemlich gut und war nicht so kompliziert wie ich erwartet hatte. Ergebnisse werden schnell generiert, haben eine physikalische Bedeutung und sehen fantastisch aus.

In diesem Artikel werde ich über meine einfache C ++ - Implementierung eines partikelbasierten hydraulischen Erosionssystems mit quadratischem Netz sprechen. Ich werde alle physikalischen Rechtfertigungen erklären, die der Implementierung zugrunde liegen, und über Mathematik sprechen. Der Code ist extrem einfach (nur etwa 20 Zeilen für die Erosionsmathematik) und schnell zu implementieren. Ich empfehle ihn daher jedem, der den Realismus seines Geländes verbessern möchte.

Die Ergebnisse werden mit einer abgespeckten Version meiner Homebrew OpenGl Engine gerendert , die ich geändert habe, um ein 2D-Array von Punkten als Höhenkarte zu rendern. Eine abgespeckte Version der Engine ist viel einfacher zu verstehen, wenn Sie OpenGL in C ++ lernen möchten.

Implementierung


Inspiriert von vielen Quellen hydraulischer Erosion entschied ich, dass es am logischsten wäre, Partikel zu verwenden.

Die Erosion auf Partikelbasis ist sehr einfach:

  • Wir erzeugen ein Teilchen an einem zufälligen Punkt auf der Oberfläche.
  • Es bewegt / gleitet entlang der Oberfläche unter Verwendung der klassischen Standardmechanik (wir werden weiter unten darĂĽber sprechen).
  • Wir fĂĽhren den Transfer von Materie / Sediment zwischen der Oberfläche und dem Partikel durch (dies wird auch unten erklärt).
  • Wir verdampfen einen Teil des Partikels
  • Wenn sich das Partikel auĂźerhalb der Karte befindet oder zu klein ist, zerstören Sie es
  • Wiederholen Sie den Vorgang mit der gewĂĽnschten Partikelanzahl.

Hinweis: Systemparameter werden im entsprechenden Abschnitt erläutert. Am wichtigsten ist der Zeitschrittfaktor dt, der alle Parameter proportional skaliert. Dies ermöglicht es uns, die Frequenz der Simulation zu erhöhen (auf Kosten der Erhöhung des Rauschens), ohne den relativen Maßstab der Parameter zu ändern. Dies ist im folgenden Code zu sehen.

Partikel


Dazu habe ich eine einfache Partikelstruktur erstellt, die alle Eigenschaften enthält, die ich benötige:

struct Particle{
  //Construct particle at position _pos
  Particle(glm::vec2 _pos){ pos = _pos; }

  glm::vec2 pos;
  glm::vec2 speed = glm::vec2(0.0); //Initialize to 0

  float volume = 1.0;   //Total particle volume
  float sediment = 0.0; //Fraction of volume that is sediment!
};

Hinweis: Falls Sie mit der GLM-Bibliothek nicht vertraut sind: Ich verwende sie, um Vektoroperationen auszufĂĽhren.

Ein Partikel hat eine Position und Geschwindigkeit, die bestimmen, wie es sich bewegt. Darüber hinaus hat es ein Volumen und einen Anteil, der bestimmt, wie viel des Volumens Sedimentgesteine ​​sind.

Uhrwerk: klassische Mechanik


Die Bewegung der Partikel auf der Oberfläche wird mithilfe der klassischen Mechanik simuliert . Kurz gesagt, die Position x des Teilchens wird durch die Geschwindigkeit v geändert, die durch die Beschleunigung a geändert wird .



Hinweis: Fettgedruckte Buchstaben zeigen an, dass der Wert ein Vektor ist.

Wir wissen auch, dass Kraft gleich Masse mal Beschleunigung ist:


Das Teilchen erfährt eine durch die Schwerkraft verursachte Abwärtsbeschleunigung, befindet sich jedoch an der Oberfläche, was eine Abwärtsbeschleunigung unmöglich macht. Stattdessen wird das Teilchen stattdessen einer Kraft F ausgesetzt, die entlang der Oberfläche gerichtet und proportional zur Normalen zur Oberfläche ist.

Daher können wir sagen, dass die Beschleunigung a proportional zum Normalenvektor der Oberfläche n geteilt durch die Masse des Teilchens ist.


Dabei ist k die Proportionalitätskonstante und m die Partikelmasse. Wenn die Masse gleich dem Volumen multipliziert mit der Dichte ist, erhalten wir das vollständige System der Teilchenbewegung unter Verwendung der klassischen Mechanik:

//... particle "drop" was spawned above at random position

glm::ivec2 ipos = drop.pos; //Floored Droplet Initial Position
glm::vec3 n = surfaceNormal(ipos.x, ipos.y);  //Surface Normal

//Accelerate particle using classical mechanics
drop.speed += dt*glm::vec2(n.x, n.z)/(drop.volume*density);
drop.pos   += dt*drop.speed;
drop.speed *= (1.0-dt*friction);  //Friction Factor

//...

Hinweis: Die Geschwindigkeit nach der Partikelbewegung wird durch den Reibungsvektor verringert. Beachten Sie, dass der Zeitschrittfaktor hier enthalten ist. Teilchen haben ihre eigene Trägheit, proportional zu ihrer Dichte, da wir ihre Bewegung mithilfe der Beschleunigung simulieren.

Sedimentationsprozess: StoffĂĽbergang


Der Prozess der Bildung von Sedimentgesteinen erfolgt physikalisch als Ăśbertragung von Sedimentgesteinen von der Erde auf das Partikel und umgekehrt an dem Punkt, an dem sich das Partikel befindet (" StoffĂĽbergang ").

In der chemischen Technologie wird der Stoffübergang (d. H. Änderung der Masse / des Sediments über die Zeit) zwischen zwei Phasen (in diesem Fall der Erdoberfläche und einem Tropfen) üblicherweise unter Verwendung von Stoffübergangskoeffizienten beschrieben .

Der StoffĂĽbergang ist proportional zur Differenz zwischen Konzentration c und Gleichgewichtskonzentration c_eq:


Dabei ist k die Proportionalitätskonstante (Stoffübergangskoeffizient). Dieser Unterschied zwischen Gleichgewicht und tatsächlicher Konzentration wird oft als "treibende Kraft" bezeichnet.

Ist die Gleichgewichtskonzentration höher als der Strom, absorbiert das Partikel Sedimentgesteine. Wenn niedriger, dann verliert sie. Wenn sie gleich sind, werden keine Änderungen vorgenommen.

Der StoffĂĽbergangskoeffizient kann auf verschiedene Arten interpretiert werden:

  • Als Häufigkeit des Ăśbergangs zwischen Phasen (hier ist es "Abscheidungsrate")
  • Die Geschwindigkeit, mit der die Tröpfchenkonzentration zur Gleichgewichtskonzentration neigt

Das System ist vollständig von der Bestimmung der Gleichgewichtskonzentration abhängig. Abhängig von der Definition zeigt das System eine unterschiedliche Dynamik der Sedimentgesteinsablagerung. In meiner Implementierung ist die Gleichgewichtskonzentration höher, wenn wir uns nach unten bewegen und wenn wir uns schneller bewegen, und sie ist auch proportional zum Volumen des Partikels:

//...

//Compute Equilibrium Sediment Content
float c_eq = drop.volume*glm::length(drop.speed)*(heightmap[ipos.x][ipos.y]-heightmap[(int)drop.pos.x][(int)drop.pos.y]);

if(c_eq < 0.0) c_eq = 0.0;

//Compute Capacity Difference ("Driving Force")
float cdiff = c_eq - drop.sediment;

//Perform the Mass Transfer!
drop.sediment += dt*depositionRate*cdiff;
heightmap[ipos.x][ipos.y] -= dt*drop.volume*depositionRate*cdiff;

//...

Hinweis: Die Änderung der Konzentration innerhalb des Partikels wird vollständig durch die Stoffübergangsgleichung beschrieben. Die Änderung der Höhenkarte wird zusätzlich mit dem Partikelvolumen multipliziert, da wir sie proportional nicht zur Konzentration, sondern zur Masse ändern (Konzentration wird mit Volumen multipliziert).

Andere Aspekte


Die Höhenkarte wird durch ein mehrschichtiges Perlin-Rauschen mit einem zufälligen Startwert initialisiert.

Am Ende jedes Zeitschritts verliert das Partikel entsprechend der Verdampfungsrate eine kleine Masse:

//...

drop.volume *= (1.0-dt*evapRate);

//...

Dieser Vorgang wird für Tausende von Partikeln wiederholt, die an zufälligen Orten erstellt und separat simuliert wurden (in meinem Fall werden die Berechnungen nacheinander in der CPU durchgeführt).

Meine Standardimplementierung enthält einen guten Satz von Parametern. Nach 200.000 Partikeln sieht die Erosion sehr gut aus.

Ergebnisse


Der fertige Code für den Erosionsprozess besteht aus ungefähr 20 Zeilen ohne Kommentar.

Hier ist ein Vorher-Nachher-Vergleich für 10 Proben. Die Simulation erzeugt sehr schöne Grate auf Höhen, Sedimentgesteine ​​lagern sich an den Seiten einiger Grate ab, was zur Schaffung schöner Plateaus führt.


Eine Auswahl von zehn Vorher-Nachher-Vergleichen. Der Shader, den ich geschrieben habe, erhält am Eingang zwei Farben. Ich habe auch Shadow Mapping, entfernungsabhängigen Nebel und Phong Shading implementiert.

Hier sind zehn weitere (verschiedene) Beispiele nur mit den Ergebnissen:


Zehn weitere Probenergebnisse. Sie unterscheiden sich von den vorherigen, auch wenn einige Formationen ähnlich aussehen.

Indem Sie den Startwert während der Karteninitialisierung unterschiedlich auswählen, können Sie kontrollierte unterschiedliche Ergebnisse erstellen.

Simulationszeit


Die Simulationszeit steht in direktem Zusammenhang mit der Lebensdauer der Partikel und der Anzahl der simulierten Partikel.

Verschiedene Faktoren beeinflussen die Lebensdauer von Partikeln und können für jedes Partikel stark variieren, da sie an zufälligen Orten erzeugt werden:

  • Reibung und Trägheit: Teilchengeschwindigkeit
  • Gittergröße: Die Wahrscheinlichkeit, dass Partikel aus der Karte fallen.
  • Verdampfungsrate: Partikelauslöschungsrate

Mit den Standardparametern benötigt eine Simulation eines einzelnen Partikels 10 bis 100 Millisekunden, was für eine Simulation von 200.000 Partikeln 10 bis 20 Sekunden ergibt.

Sie können den Erosionsgrad erhöhen und die Lebensdauer der Partikel verkürzen, indem Sie den Zeitschritt erhöhen, ohne die Anzahl der simulierten Partikel zu ändern. Dies kann die Simulation „verrauscht“ machen. Wenn der Zeitschritt zu groß ist, besteht die Möglichkeit, dass die Simulation fehlschlägt.

Die Simulation selbst in der Engine verursacht zusätzliche Kosten für die Neuerstellung des Oberflächennetzes für das Rendern.

Sie können den Code optimieren, indem Sie die Kosten für die Neuerstellung des Netzes senken oder die Geschwindigkeit eines Zeitschritts für ein einzelnes Partikel erhöhen.

Arbeit fĂĽr die Zukunft


In Kürze werde ich diesen Code in mein Territory-Projekt einfügen, um eine Höhenkarte im Höhengenerator zu erstellen.

Ich habe bereits ein Framework für die vereinfachte Simulation der Fluiddynamik in Klimasystemen geschrieben, das jedoch noch nicht veröffentlicht werden kann. Dieses System erhält Klimamuster von einer Höhenkarte. Dies wird das Thema für einen zukünftigen Beitrag über physikalisch genaue prozedurale Klimasysteme sein!

Mit einem simulierten Klimasystem wird es in Zukunft möglich sein, Partikel aus der Verteilung (z. B. an Orten, an denen es regnet) zu entnehmen, anstatt sie gleichmäßig auf der Karte zu verteilen. Sie werden ideal mit der Geologie und dem Klima des Reliefs kombiniert.

Darüber hinaus können Steine ​​nach dem Hinzufügen verschiedener Arten von geologischen Formationen einen unterschiedlichen Löslichkeitsgrad aufweisen. Dies wirkt sich direkt auf den Stoffübergangskoeffizienten (Ablagerungsrate) aus und führt zu einer anderen Erosionsrate.

Dieses System kann keinen realen Flüssigkeits- und Flussfluss simulieren, kann jedoch möglicherweise an solche Aufgaben angepasst werden. Ich werde darüber nachdenken und vielleicht werde ich in Zukunft eine Fortsetzung des Beitrags veröffentlichen.

All Articles