Normalerweise stützen sich die Sandkörner gegenseitig ab: Das Korngerüst trägt die Last. Steigt aber der Druck des Wassers in den Poren stark an, werden die Körner auseinandergedrückt. Dann überträgt das Gerüst keine Kräfte mehr, und der Boden fließt. Auslöser können Erdbeben sein, schnelle Wasserspiegeländerungen oder das Wiederbenetzen eines teilgesättigten Sandes.
Wenn Sand nass wird
Warum rutscht ein Sandhang plötzlich ab, wenn das Grundwasser steigt oder ein Schiff vorbeifährt? Wir simulieren das am Computer – vom einzelnen Korn bis zur ganzen Bodensäule.
Ein kurzer Rundgang in sechs Stationen, etwa 25 Minuten. Mit Videos aus unseren Simulationen, Reglern zum Ausprobieren und einem Selbsttest am Ende.
Eine Reise über drei Größenordnungen
Wähle einen Einstieg. Jede Skala erzählt einen Teil der Geschichte.
Die Böschung
Ufer und Kippen können sich verflüssigen. Worum geht es eigentlich?
Kapitel 1: Phänomen →Die Bodensäule
Ist der Boden hier fest oder flüssig? Ein Phasenfeld entscheidet.
Kapitel 4: Kontinuumsmodell →Der Porenraum
Wasserbrücken halten Körner zusammen. Wir simulieren jedes einzelne Korn.
Kapitel 3: Kornmodell →Wenn Boden zu fließen beginnt
Locker gelagerter, wassergesättigter Sand kann sich plötzlich wie eine Flüssigkeit verhalten. Das nennt man Verflüssigung.
≈ 4 MinutenWusstest du?
Eine Sandburg hält nur mit feuchtem Sand. Zwischen den Körnern sitzen winzige Wasserbrücken, sogenannte Menisken. Ihre Oberflächenspannung zieht die Körner zusammen: die kapillare Kohäsion. Trockener Sand rieselt, und unter Wasser verschwinden die Brücken. Genau dieser Effekt spielt in Bergbaugebieten eine große Rolle.
Sackung beim Grundwasserwiederanstieg
In ehemaligen Tagebaugebieten, etwa in der Lausitz, wurden riesige Mengen Sand locker verkippt. Nach dem Ende des Bergbaus steigt das Grundwasser wieder an.
Dabei füllen sich die Poren, die Menisken verschwinden und mit ihnen die kapillare Kohäsion. Das lockere Gerüst kann in sich zusammensacken (englisch wetting collapse). Ist der Sand danach wassergesättigt, kann er sich verflüssigen, es kommt zum Setzungsfließen.
Probiere es aus: Lass das Grundwasser steigen und beobachte die Wasserbrücken.
Grundwasser steigt
schematischVerflüssigung durch schiffsinduzierten Absunk
Wähle eine Markierung im Bild.
Die Bundesanstalt für Wasserbau (BAW) untersucht diesen Vorgang im Labor an Bodensäulen. Mit genau diesen Versuchen prüfen wir unser Kontinuumsmodell, mehr dazu in Kapitel 5.
Zwei Skalen, ein Werkzeug
Wie verteilen sich Wasser und Luft zwischen den Körnern?
Jedes Korn wird aufgelöst. Wir sehen Menisken, eingeschlossene Luft und Körner, die von Kapillarkräften verschoben werden.
Zum Kornmodell (Kapitel 3) →Wo ist der Boden fest, wo flüssig?
Keine einzelnen Körner mehr, sondern gemittelte Größen: Feststoffanteil, Porenwasserdruck und ein Phasenfeld für „sedimentiert“ oder „verflüssigt“.
Zum Kontinuumsmodell (Kapitel 4) →Beide Modelle verwenden dieselbe mathematische Idee: das Phasenfeld. Darum geht es im nächsten Kapitel.
Das Werkzeug: Phasenfelder
Statt eine Grenzfläche als scharfe Linie zu verfolgen, beschreiben wir sie durch ein glattes Feld \phi, das stetig zwischen 0 und 1 wechselt.
≈ 4 MinutenEine diffuse Grenzfläche
ausprobierenFahre mit der Maus über das Profil links oder tippe darauf, um eine Stelle x auszuwählen.
Das Double-Well-Potential W(\phi) rechts (auf Deutsch: Doppelmulde) hat zwei Täler bei \phi=0 und \phi=1: Reine Phasen kosten nichts, Mischzustände kosten Energie. Je kleiner \varepsilon, desto schmaler die Übergangszone. Für \varepsilon\to 0 nähert sich das Profil der scharfen Grenzfläche (gestrichelt). Warum der Übergang gerade etwa \varepsilon breit wird, zeigt der nächste Abschnitt.
Die Idee dahinter: eine Energie
Das Phasenfeld strebt den Zustand mit der kleinsten Energie an. Diese Energie hat zwei Anteile, die gegeneinander arbeiten. Klicke auf einen Anteil:
Der Kompromiss: Wie breit wird der Übergang?
ausprobierenDie Cahn–Hilliard-Gleichung
Wie findet \phi dieses Minimum? Die Cahn–Hilliard-Gleichung beschreibt, wie sich \phi in der Zeit verändert. Klicke auf einen Begriff:
\mu ausgerechnet Numerik
Bildet man die Variationsableitung, wird aus der Mischungsenergie die Ableitung W'(\phi) und aus der Gradientenenergie, nach partieller Integration, der Laplace-Term:
In dieser Form steht \mu in den Gleichungen von Kapitel 3 und, mit anders skaliertem Double-Well-Potential W'(\phi)=8\phi^3-12\phi^2+4\phi, im Code von Kapitel 4. Weil W' nichtlinear ist, wird es dort um den alten Zeitschritt linearisiert.
Warum der Aufwand?
Topologie geschenkt
Menisken verschmelzen, Luftblasen schnüren sich ab: Das passiert im Phasenfeld von selbst, ohne Sonderbehandlung im Code.
Eine Gleichung überall
Die Grenzfläche muss nicht explizit verfolgt werden. Sie ist einfach dort, wo \phi\approx\tfrac12 ist.
Physik ist eingebaut
Die Gleichung folgt aus einer Energie. Oberflächenspannung und Benetzung ergeben sich daraus ganz natürlich.
Was bedeutet \phi in den beiden Modellen? Überblick
Kornskala: \phi=1 ist Wasser, \phi=0 ist Luft. Das Phasenfeld lebt nur im Porenraum. Die Körner werden nicht über \phi beschrieben, sondern scharf durch ein mitbewegtes Gitter aufgelöst.
Kontinuumsskala: \phi=1 ist sedimentierter Boden (das Korngerüst trägt), \phi=0 ist verflüssigte Suspension. Hier ändert sich \phi zusätzlich durch Quellterme, denn der Boden wechselt seinen Zustand. Die Gesamtmenge von \phi bleibt dort deshalb nicht erhalten (Kapitel 4).
Mobilität: ein Detail mit großer Wirkung Numerik
Mit konstanter Mobilität tauschen Menisken unterschiedlicher Krümmung Wasser durch Diffusion aus, quer durch Luft und Wasser. Auf der Kornskala ist das ein numerisches Artefakt. Im Kornmodell verwenden wir deshalb eine entartete Mobilität, die nur in der Grenzschicht wirkt,
und schneiden \phi nach jedem Zeitschritt auf [0,1] ab. Lange Rechnungen sind ohne dieses Abschneiden gescheitert.
Das Kornmodell: jedes Korn zählt
Wir berechnen, wie Wasser und Luft durch den Porenraum strömen. Die Körner werden dabei mitgeschoben und mitgezogen.
≈ 6 MinutenStokes-Gleichungen
Im Porenraum ist alles so klein und langsam, dass Trägheit keine Rolle spielt. Es zählen nur Reibung, Druck und Oberflächenkräfte.
Cahn–Hilliard
Ein Phasenfeld \phi trennt Wasser von Luft, wie in Kapitel 2.
Mitbewegtes Gitter (ALE)
Die Kornoberflächen sind scharf im Rechengitter aufgelöst. Das Gitter bewegt sich mit den Körnern.
Die Gleichungen
Ein Zeitschritt im mitbewegten Gitter
Schritt für SchrittBenetzung: der Kontaktwinkel
ausprobierenYoung-Gesetz: Der Winkel ergibt sich aus dem Gleichgewicht der Oberflächenenergien von trockenem (\sigma_0) und benetztem Korn (\sigma_1) und der Oberflächenspannung des Wassers (\sigma).
Im Modell steckt das in einer einzigen Randbedingung, \tilde\sigma\varepsilon\,\mathbf n\cdot\nabla\phi=-\sigma_s'(\phi). Jede Höhenlinie von \phi trifft das Korn dann automatisch im richtigen Winkel.
Stimmt das auch? Ein Test mit bekannter Lösung

Je weniger Wasser in der Zelle ist, desto stärker sind die Menisken gekrümmt und desto größer ist der Kapillardruck. Die Simulation liegt auf der analytischen Kurve.
Für Neugierige
Starre Körner ohne Zusatzaufwand Numerik
Ein Korn ist im Modell eine sehr zähe Flüssigkeit, \eta_g = 10^4 bis 10^6\,\eta_w. Für \eta_g\to\infty verschwindet der symmetrische Geschwindigkeitsgradient im Korn, und \mathbf v wird zur Starrkörperbewegung. Verschiebung und Drehung folgen aus dem Kräftegleichgewicht an der Oberfläche. Es gibt also keine zusätzlichen Unbekannten. Der Preis ist eine schlecht konditionierte Matrix, deshalb verwenden wir einen direkten Löser.
Kontakt und Reibung zwischen Körnern Numerik
Ein angepasstes Gitter braucht immer einen kleinen Spalt zwischen den Körnern. Für jeden Punkt der Oberfläche wird der Abstand d zum nächsten Nachbarkorn bestimmt (k-d-Baum). Unterhalb von d_0 wirkt eine Abstoßung
Explizit ausgewertet würde sie den Zeitschritt stark beschränken. Deshalb wird der neue Abstand vorhergesagt, d^n\approx d^{n-1}+\tau\,\mathbf v^n\cdot\mathbf n, und linearisiert:
Der zweite Term dämpft die Annäherung und landet in der Matrix. Er wirkt nur, solange sich die Oberflächen nähern, damit Körner nicht aneinander kleben. Im Kontaktbereich bremst außerdem eine viskose Reibung \mathbf f_{\rm fric}=-\gamma\,(\mathbf I-\mathbf n\otimes\mathbf n)\,\mathbf v das Gleiten und Rollen.
Diskretisierung in Raum und Zeit Numerik
- Finite Elemente auf Dreiecken, die an die Kornoberflächen angepasst sind. Jede Kornoberfläche ist ein geschlossener Polygonzug aus Gitterkanten.
- Geschwindigkeit stetig, stückweise quadratisch; Druck stückweise linear mit getrennten Freiheitsgraden auf beiden Seiten der Kornoberfläche, weil er dort springt.
- \phi und \mu stückweise quadratisch, nur im Porenraum; Krümmung und Oberflächenkräfte linear auf den Kornoberflächen.
- Pro Zeitschritt ein lineares Gleichungssystem für alle Unbekannten. Nichtlinearitäten werden um den alten Zeitschritt linearisiert, z. B. \nabla\phi^{n}\otimes\nabla\phi^{n-1}.
- Oberflächenkräfte werden auf der neuen Oberfläche \mathbf x^n=\mathbf x^{n-1}+\tau\mathbf v^n ausgewertet. Das erlaubt deutlich größere Zeitschritte.
- Das Gitter wird entlang der Wasser-Luft-Grenzfläche adaptiv verfeinert.
Blick in den Code C++ · AMDiS/DUNE
Die Hauptschleife in src/ale_model.cc, stark gekürzt. Die Funktionsnamen sind die echten:
while (time < endTime) { // nach jeder Neuvernetzung von vorn NavierStokesCahnHilliard nschProb(...); // Stokes–Cahn–Hilliard + Körner Laplace laplaceProb(...); // harmonische Fortsetzung -> Gitter while (adaptInfo.time() < adaptInfo.endTime()) { timeHandler.adaptTimestepIfNecessary(); timeHandler.iteration(kappaVecOld, y0); // ① ein lineares System lösen timeHandler.moveGrid(coords, coordinates); // ② Knoten um tau*w verschieben writeFiles(...); // Ausgabe für ParaView if (remeshing(nschProb, gridPtr, ...)) // ③ Winkel zu klein? ④ neu vernetzen break; }}
Weitere Stellen: src/TimeHandler.hpp (Zeitschritt), src/NavierStokesCahnHilliard.impl.hpp (Operatoren), amdis/ALEModelAmdis2/modelExtensions/ExtendedGrains.hpp (Kontaktkräfte, Imbibition/Drainage), amdis/ALEModelAmdis2/remeshing/remeshingLoop.hh (Neuvernetzung).
Das Kontinuumsmodell: fest oder flüssig?
Auf der Skala einer Bodensäule interessiert uns nicht mehr jedes Korn. Die Frage ist: Wo trägt der Boden noch, und wo fließt er?
≈ 5 MinutenEin Phasenfeld \phi unterscheidet sedimentierten Boden (\phi=1), in dem das Korngerüst trägt, von verflüssigter Suspension (\phi=0). Je nach Wert von \phi verhält sich das Material wie ein Festkörper oder wie eine Flüssigkeit. Ausführlich beschrieben ist das Modell in Keese et al. (2026), siehe Literatur in Kapitel 6.
Neu gegenüber der klassischen Cahn–Hilliard-Gleichung sind zwei Quellterme, die von Indikatorfunktionen ein- und ausgeschaltet werden:
- \chi_{\rm liq}=1, sobald das Korngerüst keine effektiven Druckspannungen mehr überträgt: Der Boden verflüssigt.
- \chi_{\rm sed}=1, sobald der Feststoffanteil c eine Schwelle überschreitet und sich das Gerüst verdichtet: Der Boden resedimentiert.
Klicke im Ablaufdiagramm auf die Kästen.
Probier’s aus: eine Säule verflüssigen
rechnet liveHier wird die Phasenfeld-Gleichung des Modells live gelöst, eindimensional über die Tiefe. Nur die Indikatoren gibst du selbst vor. In der echten Simulation folgen sie aus Spannung und Feststoffanteil.
Von der Gleichung zum Code
Zum Schluss ein kleiner Vorgeschmack auf die Programmierung. Gerechnet wird mit AMDiS, einer C++-Bibliothek für Finite Elemente, die auf DUNE aufbaut. Jede Zeile unten hängt einen Term der schwachen Form an das Gleichungssystem an. Ausführen musst du hier nichts. Klicke auf einen Term, dann leuchten die zugehörigen Zeilen auf. Quelle: src_liquefaction/operatorFunctions.hpp.
auto opPhiTime = makeOperator(tag::test_trial{}, invTau);prob.addMatrixOperator(opPhiTime, _phi, _phi); auto opPhiOld = makeOperator(tag::test{}, invTau*phi);prob.addVectorOperator(opPhiOld, _phi); auto opPhiSource1 = makeOperator(tag::test_trial{}, k*(chi_c+chi_sigma_el));prob.addMatrixOperator(opPhiSource1, _phi, _phi); auto opPhiSource2 = makeOperator(tag::test{}, k*chi_c);prob.addVectorOperator(opPhiSource2, _phi); auto opPhiDiffusion = makeOperator(tag::partialtest_partialtrial{}, M);prob.addMatrixOperator(opPhiDiffusion, _phi, _mu); // chemical potentialauto opMu = makeOperator(tag::test_trial{}, 1.0);prob.addMatrixOperator(opMu, _mu, _mu); auto opMuLin = makeOperator(tag::test_trial{}, -(24.0*pow<2>(phi) - 24.0*phi + 4.0)/eps);prob.addMatrixOperator(opMuLin, _mu, _phi); auto opMuLinExpl = makeOperator(tag::test{}, (-16.0*pow<3>(phi) + 12.0*pow<2>(phi))/eps);prob.addVectorOperator(opMuLinExpl, _mu); auto opMuDiff = makeOperator(tag::partialtest_partialtrial{}, -eps);prob.addMatrixOperator(opMuDiff, _mu, _phi);
Die Unbekannten des Modells Überblick
- c: Volumenanteil der Feststoffphase
- p: Porenwasserdruck
- \mathbf v_s,\ \mathbf v_f: Geschwindigkeit von Feststoff und Fluid
- \phi,\ \mu: Phasenfeld und chemisches Potential
- \mathbf B: linker Cauchy-Green-Tensor (Verformung des Gerüsts), dazu elastische und viskose Spannungen
Wichtige Parameter: Sedimentationsschwelle c_{\rm sed}, Relaxationsrate K, Mobilität M, Grenzflächendicke \varepsilon.
Numerik des 1D-Säulenmodells Numerik
- Eindimensionales Gitter über die Tiefe (
Dune::OneDGrid), quadratische Lagrange-Elemente für alle neun Felder. - Implizites Euler-Verfahren; pro Zeitschritt ein gekoppeltes lineares System für Phasenfeld, Geschwindigkeiten und Mechanik.
- Adaptive Verfeinerung abhängig vom Abstand zur Grenzfläche, in drei Stufen: Grenzfläche, Übergang, Rest.
- Für Tests lässt sich die Rückwirkung des Phasenfelds auf Verformung und Spannungen abschalten (Parameter
couple_mechanics). Dann hängen Verformung und Spannungen nicht von \phi ab.
src_liquefaction/liquefaction-1d.cpp · src_liquefaction/operatorFunctions.hpp
Was die Modelle zeigen
Wie gut passt das Kontinuumsmodell zum Laborversuch? Und was passiert im Porenraum, wenn Wasser eindringt oder abfließt?
≈ 5 MinutenDer Säulenversuch der BAW
In der Wechseldurchströmungsanlage der Bundesanstalt für Wasserbau steht eine Bodensäule unter Wasser. Der Wasserdruck oben wird innerhalb von 5 Sekunden um 6 kPa abgesenkt, das entspricht etwa 0,6 m Wassersäule. So wird nachgebildet, wie der Wasserspiegel am Ufer absinkt, wenn ein Schiff vorbeifährt (schiffsinduzierter Absunk, siehe Kapitel 1). Sensoren in verschiedenen Tiefen messen den Porenwasserüberdruck \Delta p.
Erreicht \Delta p die gestrichelte Verflüssigungsgrenze (die anfängliche effektive Vertikalspannung), tragen die Körner nichts mehr: Der Boden ist dort verflüssigt.
Porenwasserüberdruck über die Tiefe
ZeitschieberDaten: Keese et al. (2026), Abb. 4 und 7. Die Messung ist ausführlich in Keese et al. (2025) beschrieben.
Im Porenraum: Simulationen des Kornmodells



Warum das wichtig ist
Kapillardruck, Sättigung und die Form der Wasserverteilung sind im Feld kaum messbar. Kornskalige Simulationen liefern sie, und zwar während sich der Porenraum durch Kornumlagerungen ständig verändert. Daraus sollen Beziehungen für das Kontinuumsmodell entstehen: Hochskalierung.
Was ist hängen geblieben?
Sieben kurze Fragen. Du siehst sofort, ob es stimmt, und warum.
Wie geht es weiter?
In die dritte Dimension
Bisher rechnen wir in 2D. Dort schneidet ein Wasserpfad die Luft darunter ab, Luft wird also leichter eingeschlossen als in echtem Sand. Die 3D-Erweiterung läuft.
Hysterese verstehen
Beim Füllen und Entleeren sieht die Wasserverteilung bei gleicher Sättigung verschieden aus. Minkowski-Funktionale (Fläche, Grenzflächenlänge, Euler-Charakteristik) machen das messbar.
Fließrutschungen
Mit reibungsbasierten Stoffmodellen soll das Kontinuumsmodell auch großräumige Fließrutschungen an Böschungen abbilden.
Zum Weiterlesen Literatur
- M. Mokbel, D. Mokbel, S. Liese, C. A. Weber, S. Aland: A simulation method for the wetting dynamics of liquid droplets on deformable membranes. SIAM J. Sci. Comput. 46 (2024), B806–B829.
- H. Keese, J. Rothschink, O. Stelzer, T. Nagel: Transient hydro-mechanical analyses of column experiments on wave-induced soil liquefaction. Comput. Geotech. 185 (2025), 107321. doi
- H. Keese, L. Zieger, O. Stelzer, S. Aland, T. Nagel: Usage of the phase field method for describing wave-induced liquefaction. In: A. Pistrol, H. F. Schweiger (Hrsg.): Proceedings of the 21st ICSMGE. ÖGG, Austrian Society for Geomechanics, 2026. doi
- J. Kayser, F. Karl, D. Schürenkamp, N. Schwab, H. Oumeraci: A test apparatus for alternating flow in geotechnical engineering. Geotech. Test. J. 39 (2016), 865–878. doi
- S. Aland, D. Mokbel: A unified numerical model for wetting of soft substrates. Int. J. Numer. Methods Eng. 122 (2021), 903–918.
- H. Tanaka, T. Araki: Simulation method of colloidal suspensions with hydrodynamic interactions: fluid particle dynamics. Phys. Rev. Lett. 85 (2000), 1338–1341.
Bild- und Videonachweis Quellen
- Simulationsvideos und -bilder (Kapitel 0, 3 und 5): eigene Simulationen, Institut für Numerische Mathematik und Optimierung, TU Bergakademie Freiberg.
- Porenwasserüberdruck im Säulenversuch (Kapitel 5): Messdaten der BAW, Referenzlösung und Phasenfeldsimulation aus Keese et al. (2026), Abb. 4 und 7. Die Messung ist ausführlich in Keese et al. (2025) beschrieben.
- Förderlogo: Europäische Union / Freistaat Sachsen.
- Alle übrigen Grafiken und Schemazeichnungen: eigene Darstellung.
Kofinanziert von der EU und mitfinanziert durch Steuermittel auf der Grundlage des vom Sächsischen Landtag beschlossenen Haushaltes.