[Python] Flüssigkeitssimulation: Inkompressible Navier-Stokes-Gleichung

Einführung

Ich möchte das Wissen zusammenfassen, das für die Erstellung eines numerischen Fluidanalysecodes für Flüssigkeiten erforderlich ist, und gleichzeitig die numerische Fluiddynamik (CFD) untersuchen. Es scheint, dass es viele Fehler gibt. Bitte kontaktieren Sie uns, wenn Sie diese finden.

Zielgruppe

Serie

Grober Inhalt dieses Artikels

Dieser Artikel befasst sich mit einphasigen numerischen Simulationen wie nur Gas oder nur Flüssigkeit. Nach einer kurzen Erläuterung der maßgeblichen Gleichung des Fluids wird es konfiguriert, um den numerischen Berechnungsalgorithmus des inkompressiblen Fluids zu implementieren, der für die Simulation des Fluids erforderlich ist. Schließlich simulieren wir den Hohlraumfluss (Gefühl des Schleuderns der Flüssigkeit im Bereich / Abbildung unten). In Bezug auf den Hohlraumfluss ist @ eigs 'Numerische Lösung der inkompressiblen Navier-Stokes-Gleichung 1: Einführung leicht zu verstehen und die Grafik ist auch schön.

2d-cavity.gif

Inhaltsverzeichnis

Kapitel Titel Bemerkungen
1. Flüssige Lösung
1.1. Flüssigkeitsregelgleichung
1.2. Komprimierbar und unkomprimierbar
2. Berechnungsalgorithmus für inkompressible Flüssigkeiten
2.1. MAC-Methode Marker And Cell method
2.2. Fractional Step-Methode Wird auch als Teilschrittmethode bezeichnet
2.3. EINFACHE Methode Semi-Implicit Method for Pressure-Linked Equation method
2.4. CCUP-Methode CIP-Combined Unified Procedure method
2.5. Räumlicher Unterschied In Bezug auf die variable Platzierung
3. Implementierung CCUP-Methode

1. Flüssige Lösung

1-1. Kontrollgleichung der Flüssigkeit

Eine Flüssigkeit besteht aus drei Gleichungen: einer kontinuierlichen Gleichung, einer Navier-Stokes-Gleichung und einer Energiegleichung. Jede Gleichung wird aus Massenerhaltung, Impulserhaltung und Energieeinsparung abgeleitet. Es mag verwirrender sein, aber ich werde auch ein Bild von jeder Formel hinzufügen.

a. Kontinuierlicher Ausdruck

\frac{\partial \rho}{\partial t} +  (\vec{u} \cdot \nabla) \rho  =   - \rho \nabla \cdot \vec{u}  \\
\quad or \\
\frac{\partial \rho}{\partial t} + \nabla \rho \vec{u} = 0

Es kann aus der Tatsache abgeleitet werden, dass die Massenänderung innerhalb eines bestimmten Testvolumens gleich dem Massenstrom $ \ rho \ vec {u} $ ist, der in die Oberfläche hinein und aus dieser heraus fließt, geteilt durch die Fläche. Um dies anhand eines leicht verständlichen Beispiels zu erklären, zeigt die obige Formel, dass der Betrag der Geldveränderung in der Brieftasche gleich dem Betrag ist, der durch Subtrahieren des Geldbetrags erhalten wird, der durch die Tür der Brieftasche gelangt ist.

money.png

b. Navier-Stokes-Gleichung (Fluid-Bewegungsgleichung)

\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} \\
 \quad or \\
\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \nabla \cdot {\bf T}_{\nu} + \vec{F_B} \\
, where \quad {\bf T}_{\nu} = \lambda (\nabla \cdot \vec{u}) {\bf I} + 2 \mu {\bf D} \quad :Viskoser Spannungstensor\\
{\bf D}  = \frac{1}{2} \left\{(\nabla \vec{u}) + (\nabla \vec{u})^T \right\} :Dehnungsratentensor\\
\quad \vec{F_B} \quad :Volumenkraft(Hauptsächlich Schwerkraft) \\
\quad \mu :Viskositätskoeffizient

Der Impuls der Flüssigkeit ist $ Druckgradient: \ nabla p $, $ Diffusion aufgrund der Viskosität: \ nabla \ cdot {\ bf T} _ {\ nu} $, $ Volumenkraft (Schwerkraft usw.): \ vec {F_B} $ Zeigt an, dass Sie betroffen sind. Um ein Bildbeispiel zu geben, kann es hilfreich sein, sich vorzustellen, dass sich die Geschwindigkeit des Balls aufgrund der Auswirkungen von Druckdifferenz, Viskosität und Schwerkraft ändert, wenn ein Junge einen Ball tritt.

navier_stokes.png

c. Energiegleichung

 \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L} \\
, where \quad  \dot{L} :Wärmemenge, die durch chemische Reaktionen usw. erzeugt wird.\\
\quad \kappa :Wärmeleitfähigkeit

Jedoch,

{\bf T}_{\nu} : {\bf D} = tr({\bf T}_{\nu} \cdot {\bf D})

Das Ausmaß der Änderung der inneren Energie des Fluids ist die Arbeit $ \ nabla \ cdot \ vec {u} $ aufgrund des Drucks, die Energiebilanz aufgrund der Wärmeleitung $ \ nabla \ cdot (\ kappa \ nabla T) $ und die dissipierte Energie aufgrund der Viskosität $ {\ bf Dies bedeutet, dass es von $ \ dot {L} $ wie T} _ {\ nu} \ cdot \ nabla \ vec {u} $ und chemischen Reaktionen betroffen ist. Mit dieser Formel können Sie die Temperaturänderung der Flüssigkeit berechnen. Kurz gesagt ist die Temperaturänderung des Fluids gleich der zugeführten Wärmemenge.

energy_equation.png

Zusammenfassung

Die oben gezeigte Abbildung kann wie folgt zusammengefasst werden.

flow_equations.png

Oben ist die maßgebliche Flüssigkeitsgleichung in einem nicht konservativen System dargestellt. Wenn die rechte Seite auf 0 gesetzt ist, bedeutet dies die Übertragungsgleichung von Dichte, Geschwindigkeit und innerer Energie (eine Gleichung, die die Bewegung von Substanzen bedeutet). Ein nicht konservatives System ist eine Form, in der Variablen außerhalb der Differenzierung existieren, und es gibt keine Garantie dafür, dass die Berechnungsergebnisse dem Erhaltungsgesetz entsprechen. Daher muss die Konservativität bei der Durchführung numerischer Berechnungen bestätigt werden. Das Gegenteil des nichtkonservativen Systems wird als konservatives System bezeichnet. Dies ist ein Ausdruck, bei dem alle Variablen im Differential enthalten sind.

1-2. Komprimierbar und nicht komprimierbar

Die in Abschnitt 1-1 gezeigte Gleichung ist die maßgebliche Gleichung für komprimierbare Flüssigkeiten. Ein komprimierbares Fluid ist ein Fluid, dessen Strömungsgeschwindigkeit etwa das 0,3-fache der Schallgeschwindigkeit beträgt (Machzahl ist 0,3), und ein Fluid mit einem langsameren Fluss wird allgemein als inkompressibles Fluid erkannt. Kurz gesagt, wenn sich die Flüssigkeit zu schnell bewegt, dehnt sich die Flüssigkeit selbst aus und zieht sich zusammen (Änderungen in der Dichte), wodurch sie komprimierbar wird. Wenn sie langsam oder dicht ist, dehnt sich die Flüssigkeit nicht aus oder zieht sich zusammen und wird als inkompressibel behandelt. Daher wird die Flüssigkeit grundsätzlich als nicht komprimiert behandelt.

flow.png

Die in Abschnitt 1-1 gezeigte Formel wird in eine unkomprimierte Formel umgewandelt. Alles was wir tun müssen, ist die Dichte als Konstante und nicht als Variable zu behandeln, um die Gleichung zu vereinfachen. Die Navier Stokes- und Energiegleichungen werden durch kontinuierliche Gleichungen vereinfacht. Darüber hinaus werden wir hier die Auswirkungen chemischer Reaktionen ignorieren.

a. Kontinuierlicher Ausdruck

\nabla \cdot \vec{u} = 0

b. Navier Stokes-Gleichung

\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + \vec{F_B}

c. Energiegleichung

\rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = - \nabla \cdot (\kappa \nabla T) + 2 \mu {\bf D}:{\bf D}

Die Gleichungen für kompressives und nicht komprimierbares Fluid können wie folgt zusammengefasst werden.

Grundgleichung Kompressibilität 非Kompressibilität
a.Kontinuierliche Formel \frac{\partial \rho}{\partial t} + (\vec{u} \cdot \nabla) \rho = - \rho \nabla \cdot \vec{u} $ \nabla \cdot \vec{u} = 0$
b.Navier Stokes-Gleichung $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} $ $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + \vec{F_B} $
c.Energiegleichung $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L}$ $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = - \nabla \cdot (\kappa \nabla T) + 2 \mu {\bf D}:{\bf D}$

2. Berechnungsalgorithmus für inkompressible Flüssigkeiten

Da wir Flüssigkeiten simulieren möchten, wird in diesem Kapitel beschrieben, wie inkompressible Flüssigkeiten gelöst werden. Grundsätzlich werden wir eine diskrete Gleichung unter Verwendung einer inkompressiblen Gleichung erstellen. Um Komplikationen zu vermeiden, werden wir hier die Energiegleichungen ignorieren und "explizit" mit der expliziten Methode und "implizit lösen" mit der impliziten Methode lösen. Grundsätzlich verweise ich auf Compressible Fluid Dynamics.

Typische Berechnungsalgorithmen sind in der folgenden Tabelle zusammengefasst. Die MAC-Methode (Maker And Cell), die Fractional Step-Methode und die SIMPLE-Methode sind die grundlegenden Berechnungsalgorithmen für inkompressible Flüssigkeiten. Da es in verschiedenen Sites und Büchern implementiert ist, werde ich hier nur das Prinzip erläutern. Am Ende dieses Kapitels erwähnen wir die CCUP-Methode (CIP-Combined Unified Procedure) und implementieren sie im nächsten Kapitel. Die CCUP-Methode ist eine Methode zum Lösen der Fluidgleichung mit der in den vorherigen Artikeln verwendeten CIP-Methode.

Algorithmus Translokationsbegriff Diffusionsbegriff Druckbegriff
MAC-System Explizite Methode Explizite Methode Implizite Methode
Fractional Step-Methode Explizite Methode Implizite Methode Implizite Methode
---- Implizite Methode Explizite Methode Implizite Methode
EINFACHES Methodensystem Implizite Methode Implizite Methode Implizite Methode
CCUP-Methode Explizite Methode Explizite MethodeorImplizite Methode Implizite Methode

2-1. MAC-Methode

Die MAC-Methode (Maker And Cell) ist eine grundlegende Lösung für die Flüssigkeitsberechnung, da sie verschiedene Flüsse auf stabile Weise lösen kann. Das Berechnungsverfahren ist wie folgt.

  1. Lösen Sie die Druckgleichung mit der impliziten Methode, um den Druck $ p ^ {n + 1} $ zu ermitteln.
  2. Erweitern Sie die Geschwindigkeitsgleichung im Laufe der Zeit, um die Geschwindigkeit $ \ vec {u ^ {n + 1}} $ zu ermitteln.

Die Druckgleichung und die Geschwindigkeitsgleichung werden wie folgt berechnet. Multiplizieren Sie die Navier Stokes-Gleichung mit $ \ nabla $, um die Divergenz zu nehmen und zu zerstreuen.

\nabla^2 p^{n+1} = - \nabla \cdot \left\{(\vec{u^n} \cdot \nabla )\vec{u^n} \right\} + \frac{\nabla \cdot \vec{u^n}}{\Delta t} + \nabla \cdot \left\{ \nu \nabla^2 \vec{u^n} + \frac{1}{\rho}\vec{F}_B\right\}

Es kann erhalten werden. Wenn der Translokationsterm und der Viskositätsterm der Navier-Stokes-Gleichung explizit zeitlich getrennt sind,

    \frac{\vec{u^{n+1}} - \vec{u^n}}{\Delta t} = - (\vec{u} \cdot \nabla)\vec{u^n} - \frac{1}{\rho} \nabla p^n + \nu \nabla^2 \vec{u^n} \\
    \nabla \cdot \vec{u^{n+1}} = 0

Es kann erhalten werden. Indem wir $ \ nabla \ cdot \ vec {u ^ n} $ absichtlich lösen, ohne es auf 0 zu setzen, entwickeln wir Möglichkeiten, um Fehler während der Verteilung zu unterdrücken.

Es gibt die SMAC-Methode (Simplified MAC), die HSMAC-Methode (Highly Simplified MAC) usw. als Berechnungsalgorithmus für das MAC-Methodensystem. Diese wurden vorgeschlagen, um die MAC-Methode mit hoher Geschwindigkeit zu lösen.

2-2. Fractional Step-Methode

Bei der MAC-Methode wird der viskose Term durch die explizite Methode gelöst und wird stark von der Zeitschrittbreitenbeschränkung (CFL-Bedingung) beeinflusst. Die in diesem Abschnitt beschriebene Fractional-Step-Methode löst das Problem der CFL-Bedingungen durch implizite Bewertung des Viskositätsterms. Wird auch als Teilschrittmethode bezeichnet. Als Merkmal wurden in der Formel, die sich auf die bei der MAC-Methode verwendete Geschwindigkeit bezieht, Geschwindigkeit und Druck auf der rechten Seite gemischt, aber bei der Fractional Step-Methode werden Geschwindigkeit und Druck auf der rechten Seite wie folgt vollständig getrennt.

\frac{\bar{u} - \vec{u^n}}{\Delta t} = - (\vec{u} \cdot \nabla)\vec{u^n} + \nu \nabla^2 \vec{u^n} \\
    \frac{\vec{u^{n+1}} - \bar{u}}{\Delta t} = - \frac{1}{\rho} \nabla p^{n+1} \\
    \nabla \cdot \vec{u^{n+1}} = 0
\nabla^2 p^{n+1} = \frac{\nabla \cdot \bar{u}}{\Delta t}

Als Berechnungsverfahren

  1. Ermitteln Sie die vorläufige Geschwindigkeit $ \ bar {u} $ aus der oberen Formel für die Geschwindigkeitsformel.
  2. Verwenden Sie die Druckgleichung, um den Druck $ p ^ {n + 1} $ zu ermitteln.
  3. Ermitteln Sie die Geschwindigkeit $ \ vec {u ^ {n + 1}} $ mithilfe der zweiten Formel oben in der Formel für die Geschwindigkeit.

2-3. EINFACHE Methode

Abkürzung für Semi-Implicit Method for Pressure-Linked Equation-Methode, die basierend auf der impliziten Methode formuliert wird. will mehr wissen! Grundlagen der Thermofluidanalyse 64 Kapitel 6 Thermofluidanalysemethode: 6.5.4 EINFACHE Methode wird ausführlich erläutert.

  1. Ermitteln Sie implizit die vorläufige Geschwindigkeit $ \ bar {u} $.
\bar{u} = \vec{u^n} - \Delta t \left\{(\bar{u} \cdot \nabla)\bar{u} + \nabla p - \nu \nabla^2 \bar{u}\right\}
  1. Aktualisieren Sie den Druck mithilfe der Druckgleichung.
\nabla^2 \delta p = \frac{\nabla \cdot \bar{u}}{\Delta t}\\ p^{n+1} = p + \delta p
  1. Aktualisieren Sie die Geschwindigkeit $ \ vec {u} $ mit dem berechneten Druck.
\vec{u^{n+1}} - \bar{u} = - \Delta t \nabla(p^{n+1} - p)

2-4. CCUP-Methode

Die CCUP-Methode ist eine von Yabe et al. Konstruierte Methode, die die CIP-Methode entwickelt hat, wie wenn die Burgers-Gleichung in vorheriger Artikel durch die CIP-Methode gelöst wurde. Die maßgebliche Gleichung des Fluids wird getrennt für die Translokationsphase und die Nicht-Translokationsphase berechnet. Shimizu et al. und Himeno et al. ) Das Papier ist leicht zu verstehen. Da die Reihenfolge der Lösung der Translokations- und Nicht-Translokationsbegriffe von Person zu Person unterschiedlich ist, werden wir sie in der Reihenfolge lösen, die im letztgenannten Dokument unten eingeführt wurde. Als Algorithmus

  1. Übertragen Sie in der Übertragungsphase $ Dichte: \ rho $, $ Geschwindigkeit: \ vec {u} $, $ interne Energie: e $, $ Druck: p $.
  2. Wechseln Sie in die Phase ohne Übertragung.
  3. Aktualisieren Sie Geschwindigkeit, innere Energie und Druck in der Diffusionsphase. (Dichte ist konstant)
  4. Lösen Sie die Druckgleichung implizit und aktualisieren Sie den Druck.
  5. Erhalten Sie Dichte, Geschwindigkeit und innere Energie mit neuem Druck.

Ich werde in der Reihenfolge lösen. Kurz gesagt, nach der Übertragung und Diffusion der Substanz denke ich, dass die Bewegungsgeschwindigkeit überwältigend hoch ist (die gleiche Geschwindigkeit wie die Schallgeschwindigkeit) und der Tsuji 褄 angepasst wird.

ccup_method.png

2-5. Räumlicher Unterschied

Bevor wir mit der Implementierung beginnen, können wir sie implementieren, ohne die Methode der räumlichen Differenz zu kennen, aber ich werde sie kurz erwähnen. Dies liegt daran, dass Druckvibrationen, die als Schachbrettinstabilität bezeichnet werden, auftreten können, je nachdem, wo sich die Variable im Raum befindet.

Regelmäßiges Gitter

regular_stencil.png

Versetztes Gitter

staggered_grid.png

Colocate Gitter

collocated_grid_2.png

Bei der obigen Methode besteht die Grundrichtlinie darin, das versetzte Gitter für die Berechnung zu verwenden. Bei der Implementierung dieses Artikels wird die numerische Berechnung jedoch unter Verwendung des regulären Gitters durchgeführt (~~ Ich konnte die Motivation bei Verwendung des versetzten Gitters usw. nicht verstehen. Also ~~).

3. Implementierung

Es wird durch die CCUP-Methode implementiert. Das zu lösende Problem ist der Hohlraumfluss, der sich anfühlt, als würde er der oberen Wand eine Quergeschwindigkeit verleihen, was einen Wind im Uhrzeigersinn in dem Bereich verursacht. Der Code wird auf [Github] veröffentlicht (https://github.com/kqts/flow_simulation/blob/master/2d_incompressible_flow.ipynb). Ignoriere die Energiegleichung. Da der Übertragungsterm und der Diffusionsterm im vorherigen Artikel beschrieben wurden, ist der neue Inhalt im Code der Druckgleichungsteil, aber die Berechnung des Diffusionsterms und Die Struktur ist fast gleich. Bei der Druckberechnung werden einige Übertragungsbedingungen weggelassen. Auch die Randbedingungen sind ziemlich grob

cavity_flow.png

Wenn Sie die Geschwindigkeit mit einer roten Linie und den Druck mit einer blauen Kontur ausdrücken, können Sie sehen, dass ein solches Diagramm ausgegeben wird.

2d-cavity.gif

scipy.sparse.dia_matrix

Wir verwenden scipy.sparse.dia_matrix, um die eindimensionale simultane Gleichungsmatrix $ {\ bf A} $ $ {\ bf A} \ vec {x} = \ vec {b} $ zu erstellen. Diese Matrix $ {\ bf A} $ wird als Sparse-Matrix bezeichnet. Da die meisten Elemente der Matrix aus 0 bestehen, wird sie durch scipy.sparse.dia_matrix [Diagonal Compression Storage Method (DIA)] angegeben. ](Http://www.turbare.net/transl/scipy-lecture-notes/advanced/scipy_sparse/dia_matrix.html) ist speichereffizient und so weiter. Andere Sparse-Matrix-Speicherformate werden in Python-Artikeln beschrieben, die eine Hochgeschwindigkeits-Sparse-Matrix-Berechnung ermöglichen (https://qiita.com/KQTS/items/e5500ba6e2681456e268).

Ein Beispiel ist unten gezeigt.

#Bereiten Sie ein oder mehrere Arrays gleicher Länge vor. Beachten Sie, dass Sie die gleiche Länge benötigen.
data = np.array([np.arange(10), np.arange(20,30)], np.arange(30,40))  
#Bestimmen Sie den Versatz der Diagonalmatrix anhand des Datenindex. 0:Diagonale Matrix,Positive Zahl:Oberseite,Negative Zahl:Niedriger
#Beachten Sie, dass bei einem positiven Versatz das Array von Anfang an zugeschnitten wird und bei einem negativen Versatz das Array von hinten zugeschnitten wird.
offsets = np.array([0, 2, -3])  
a_matrix = scipy.sparse.dia_matrix((data, offsets), shape=(len(data[0]), len(data[0])))  #Stellen Sie die Form der Matrix mit der Form ein.

print(a_matrix)
# <10x10 sparse matrix of type '<class 'numpy.int64'>'
# 	with 25 stored elements (3 diagonals) in DIAgonal format>

print(a_matrix.todense())
# matrix([[ 0,  0, 22,  0,  0,  0,  0,  0,  0,  0],
#         [ 0,  1,  0, 23,  0,  0,  0,  0,  0,  0],
#         [ 0,  0,  2,  0, 24,  0,  0,  0,  0,  0],
#         [30,  0,  0,  3,  0, 25,  0,  0,  0,  0],
#         [ 0, 31,  0,  0,  4,  0, 26,  0,  0,  0],
#         [ 0,  0, 32,  0,  0,  5,  0, 27,  0,  0],
#         [ 0,  0,  0, 33,  0,  0,  6,  0, 28,  0],
#         [ 0,  0,  0,  0, 34,  0,  0,  7,  0, 29],
#         [ 0,  0,  0,  0,  0, 35,  0,  0,  8,  0],
#         [ 0,  0,  0,  0,  0,  0, 36,  0,  0,  9]])

4. Zusammenfassung

Als nächstes möchte ich über die Gas-Flüssigkeits-Zweiphasensimulation (Gas und Flüssigkeit) sprechen.

Verweise

Recommended Posts

[Python] Flüssigkeitssimulation: Inkompressible Navier-Stokes-Gleichung
[Python] Flüssigkeitssimulation: Implementieren Sie die Übertragungsgleichung
[Python] Fluidsimulation: Implementieren Sie die Diffusionsgleichung
[Python] Fluidsimulation: Von linear zu nichtlinear
Ambisonics Simulation Python
Erste Nervenzellsimulation mit NEURON + Python
Ambisonics-Simulation (externes Problem) Python
Röntgenbeugungssimulation mit Python