A family of second order, linear, unconditionally stable methods for the Cahn-Hilliard-Navier-Stokes equations

Technische Zusammenfassung: Eine Familie von Verfahren zweiter Ordnung, linear und bedingungslos stabil, für die Cahn-Hilliard-Navier-Stokes-Gleichungen

Problemstellung
Die genaue numerische Simulation komplexer Mehrphasen-Fluid-Dynamiken, insbesondere jener, die durch die Cahn-Hilliard-Navier-Stokes (CHNS)-Gleichungen modelliert werden, stellt erhebliche Herausforderungen dar, die auf der mathematischen Steifigkeit des vierten Ordnung aufweisenden Cahn-Hilliard-Operators und dessen nichtlinearer Kopplung mit den inkompressiblen Navier-Stokes-Gleichungen beruhen. Standardmäßige explizite Zeitschrittverfahren unterliegen strengen Zeitschrittbeschränkungen (z. B. $\Delta t \sim O(\Delta x^4)$), um die Stabilität zu gewährleisten. Während voll implizite oder Convex-Splitting-Verfahren eine bedingungslose Energiestabilität bieten, erfordern sie oft das Lösen großer, gekoppelter nichtlinearer Systeme bei jedem Zeitschritt, was den Rechenaufwand erhöht. Es besteht ein Bedarf an Algorithmen, die eine zeitliche Genauigkeit zweiter Ordnung, bedingungslose Stabilität und lineare Lösbarkeit kombinieren, ohne die physikalischen Erhaltungseigenschaften des Systems zu opfern.

Methodik
Die Autoren schlagen eine Familie von IMEX-Verfahren (Implicit-Explicit) zweiter Ordnung, linear und bedingungslos stabil für Zwei-Phasen-Strömungen mit übereinstimmender Dichte vor. Das semidiskrete Schema integriert drei primäre Komponenten:

1. Extrapolation nichtlinearer Terme: Nichtlineare Konvektions- und Kopplungsterme werden explizit unter Verwendung extrapolierter Größen (bezeichnet als $H_{n+\theta}$) behandelt, wodurch nichtlineare Iterationen vermieden werden.
2. Hilfsvariablen-Formulierung: Der nichtlineare Freie-Energie-Term wird über eine Hilfsvariable $q$ (abgeleitet aus dem Doppelmuldenpotenzial $F(\phi)$) behandelt, was das Problem zur Aufrechterhaltung der Linearität reformuliert.
3. Zeitliche Krümmungsregularisierung: Ein Stabilisierungsmechanismus, der durch einen Parameter $\epsilon \geq 0$ gesteuert wird, wird eingeführt. Dies beinhaltet einen spezifischen Interpolationsoperator $J^\epsilon_{n+\theta}$, der die diskrete zeitliche Krümmung der Lösung einbezieht – ein Konzept, das aus vorangegangenen Arbeiten zur Navier-Stokes-Regularisierung adaptiert wurde.

Der resultierende Algorithmus (Algorithmus 2.1) erfordert bei jedem Zeitschritt lediglich das Lösen linearer Systeme mit konstanten Koeffizienten. Das Schema ist parametrisiert durch $\theta \in (1/2, 1]$ und $\epsilon \geq 0$.

Wesentliche Beiträge und theoretische Ergebnisse

• Bedingungslose Stabilität: Die Arbeit etabliert eine rigorose diskrete Energieabschätzung, die beweist, dass das vorgeschlagene Schema für $\theta \in (1/2, 1]$ und $\epsilon \geq 0$ bedingungslos langfristig stabil ist. Der Beweis nutzt symmetrische positiv definite Matrizen, um diskrete Normen zu definieren, und zeigt, dass das Schema ein diskretes Energiedissipationsgesetz ohne Zeitschrittbeschränkungen erfüllt.
• Lineare Lösbarkeit: Im Gegensatz zu voll impliziten oder Convex-Splitting-Ansätzen, die möglicherweise nichtlineare Solver erfordern, liefert diese Methode bei jedem Zeitschritt lineare Systeme, was die Komplexität der Berechnung erheblich reduziert.
• Genauigkeit zweiter Ordnung: Theoretische Analysen und numerische Experimente bestätigen, dass die Methode eine zeitliche Genauigkeit zweiter Ordnung beibehält.

Numerische Ergebnisse
Die Autoren validieren die Methode durch eine Reihe von Benchmark-Berechnungen:

• Konvergenzanalyse: Unter Verwendung der Methode der künstlich erzeugten Lösungen (Method of Manufactured Solutions, MMS) demonstriert das Schema eine annähernd zeitliche Konvergenz zweiter Ordnung für Geschwindigkeit, Druck und Phasenfeldvariablen über verschiedene Finite-Elemente-Konfigurationen ($P2-P1-P2-P2$ und $P3-P2-P3-P3$) hinweg.
• Robustheit und Regularisierung: In der Zwei-Phasen-Lid-Driven-Cavity-Problematik untersucht die Studie das Zusammenspiel von $\theta$ und $\epsilon$. Es wird beobachtet, dass das Schema für $\theta$-Werte nahe der Crank-Nicolson-Grenze (z. B. $\theta=0.51$) mit $\epsilon=0$ aufgrund schwach gedämpfter zeitlicher Oszillationen instabil werden kann. Die Einführung einer positiven Krümmungsregularisierung ($\epsilon > 0$) dämpft diese Oszillationen effektiv und stellt die Beschränktheit sowie Robustheit selbst bei größeren Zeitschritten wieder her.
• Physikalische Benchmarks:
  • Spinodale Entmischung: Die Methode erfasst die Energiedissipation und Massenerhaltung in Langzeit-Simulationen ($T=100$) präzise und zeigt eine effektive Phasentrennung und Vergröberung (Coarsening).
  • Tropfenrelaxation: Der Algorithmus löst die oberflächenspannungsgetriebene Grenzflächenbewegung erfolgreich auf, wobei er den Übergang eines quadratischen Tropfens zu einem kreisförmigen Gleichgewicht unter Erhaltung der Masse vollzieht.
  • Lid-Driven Cavity & Rayleigh-Taylor-Instabilität: Das Schema bewältigt starke Scherung, Grenzflächendeformation und dichtegetriebene Instabilitäten (einschließlich Kelvin-Helmholtz-Roll-up in Niedrigviskositätsregimen) ohne sichtbare numerische Oszillationen.

Bedeutung
Die Arbeit behauptet, dass die vorgeschlagene Familie von Methoden eine praktische und effiziente Alternative zu voll nichtlinearen Diskretisationen für CHNS-Systeme bietet. Durch die Entkopplung der Nichtlinearitäten mittels Extrapolation und Hilfsvariablen sowie die Stabilisierung der zeitlichen Diskretisation via Krümmungsregularisierung erreicht die Methode ein Gleichgewicht zwischen Genauigkeit zweiter Ordnung, bedingungsloser Stabilität und linearem Rechenaufwand. Die Autoren betonen, dass das Schema die wesentlichen physikalischen Strukturen des Modells, einschließlich Energiedissipation und Massenerhaltung, bewahrt, was es für die Simulation komplexer Grenzflächenströmungen in verschiedenen Regimen geeignet macht.