[Python] Flüssigkeitssimulation: Implementieren Sie die Übertragungsgleichung

Einführung

Ich möchte (in mehreren Artikeln) das Wissen zusammenfassen, das zum Erstellen eines Codes für die numerische Flüssigkeitsanalyse für Wasser erforderlich ist, und gleichzeitig die numerische Fluiddynamik (CFD) untersuchen, eine Studie zur Simulation von Flüssigkeiten wie Luft und Wasser.

Ich möchte es schreiben, damit auch Anfänger es leicht verstehen können. Es scheint, dass es viele Fehler gibt. Bitte kontaktieren Sie uns, wenn Sie diese finden. Ich würde mich auch freuen, wenn Sie kommentieren könnten, was schwer zu verstehen ist. Wir werden es von Zeit zu Zeit aktualisieren.

Zielgruppe

• Leute, die Python benutzen können
• Diejenigen, die an numerischen Berechnungen interessiert sind
• Diejenigen, die sich für Fluiddynamik interessieren
• Diejenigen, die grundlegende Universitätsphysik und Mathematik verstehen (über Differentialgleichungen?)

Serie

• Kapitel 1: Dieser Artikel
• Kapitel 2: [Python] Fluid Simulation: Implementieren von Diffusionsgleichungen
• Ergänzende Position in Kapitel 2: [Python] Artikel, der die Hochgeschwindigkeits-Sparse-Matrix-Berechnung ermöglicht
• Kapitel 3: [Python] Fluid Simulation: Von linear zu nicht linear
• Kapitel 4: [Python] Fluid Simulation: Inkompressible Navier-Stokes-Gleichung

Grober Inhalt dieses Artikels

Als ersten Schritt zur Behandlung der für die Wassersimulation erforderlichen Grundgleichungen von Flüssigkeit werden wir die Übertragungsgleichungen kurz zusammenfassen und implementieren.

Änderungsprotokoll

• 2020/02/19 Es gab einen Fehler im Text, deshalb habe ich ihn korrigiert. Vielen Dank an @picric_acid für den Kommentar!
• 2020/02/21 Es gab einen Fehler im Code, der behoben wurde. Vielen Dank an @fluxCapacitor für den Kommentar!

Translokationsgleichung (Formel für die Bewegung von Substanzen)

Wir werden die eindimensionale Übertragungsgleichung nach der Finite-Differenzen-Methode lösen.

Vorausgleichung $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $

Diese Gleichung bedeutet, dass sich eine Funktion u mit der Geschwindigkeit c bewegt. Kurz gesagt, es ist eine Gleichung, die ausdrückt, dass sich die physikalische Größe entlang des Geschwindigkeitsflusses c bewegt.

Zerstreue dies und löse die Gleichung. Dieses Mal werden wir es mit einer Methode namens Positive Solution lösen. Die explizite Methode ist einfach eine Methode, um zukünftige Werte nur unter Verwendung der aktuellen Zeitwerte zu finden. Neben der expliziten Methode gibt es auch eine Methode namens implizite Methode, bei der es sich um eine Berechnungsmethode handelt, die auf dem aktuellen Wert und dem zukünftigen Wert basiert.

Die explizite Methode, die dieses Mal implementiert werden soll, lautet wie folgt. Ich habe die folgenden vier nach Belieben ausgewählt.

• FTCS(Forward in Time and Central Difference in Space)
• Gegenwinddifferenzierung
• Lax-Wendroff scheme
• CIP(Constrained Interpolation Profile scheme) method

Für jeden werde ich die diskreten Ausdrücke aufschreiben. In der Formel bedeutet der Index j Ort und n Zeit. Daher bedeutet $ f_j ^ n $ den Wert der Funktion f, die sich zum Zeitpunkt $ n $ befindet, und den Ort $ x_j $.

1. FTCS(Forward in Time and Central Difference in Space) 2. $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) $

2. Upwind differencing 3. $ u_j^{n+1} = u_j^{n} - c \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j}^n - u_{j-1}^n \right) $

3. Lax-Wendroff scheme 4. $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) + \frac{c^2}{2} \left( \frac{\Delta t}{\Delta x} \right)^2 \left(u_{j+1}^n - 2 u_j^n+ u_{j-1}^n \right) $

4. CIP(Constrained Interpolation Profile scheme) method

Details werden später beschrieben. Alternativ finden Sie unter Einführung in das KVP-Recht.

Stabile Bedingungen der expliziten Methode

Ich möchte kurz auf die wichtigen Stabilitätsbedingungen bei Verwendung der expliziten Methode eingehen. Numerisch muss die Informationsausbreitungsgeschwindigkeit $ \ frac {\ Delta t} {\ Delta x} $ (die Geschwindigkeit, die durch Teilen der Gitterbreite durch die Zeit erhalten wird) größer oder gleich der physikalischen Signalgeschwindigkeit $ c $ sein. Andernfalls bewegt sich das physikalische Phänomen schneller als die numerische Berechnung, und das physikalische Phänomen kann nicht numerisch ausgedrückt werden. Basierend auf dieser Idee $ c \leq \frac{\Delta x}{\Delta t} $ Die Bedingung wird abgeleitet. Transformierte dies $ \nu = c \frac{\Delta t}{\Delta x} \leq 1 $ Wird als CFL-Bedingung bezeichnet und als numerisch stabile Bedingung verwendet. Die linke Seite wird auch als Coolan-Nummer bezeichnet. Kurz gesagt, die Zeitschrittbreite $ \ Delta t $ sollte kleiner oder gleich $ \ frac {\ Delta x} {c} $ sein. Beachten Sie, dass diese CFL-Bedingungen bei Verwendung der impliziten Methode irrelevant sind, sodass die implizite Methode den Vorteil hat, dass die Zeitinkremente im Vergleich zur expliziten Methode erhöht werden können. Die Berechnung wird jedoch entsprechend kompliziert.

Implementierung

Die verwendete Sprache ist Python. Ich werde es schreiben, ohne viele Numpy-Funktionen zu verwenden, damit die Berechnungsformel leicht zu verstehen ist.

Berechnungsziel

Berechnen Sie die Übertragung einer Rechteckwelle. Die Gitterbreite $ \ Delta x $, die physikalische Signalgeschwindigkeit $ c $ und die Zeitschrittbreite $ \ Delta t $ werden alle als Konstanten berechnet. Dieses Mal wird es gemäß der Referenz als $ \ Delta x = 1 $, $ c = 1 $, $ \ Delta t = 0,2 $ berechnet. Daher ist es auf $ CFL = 0,2 $ festgelegt und die CFL-Bedingung ist erfüllt.

# Creat square wave
Num_stencil_x = 101
x_array = np.arange(Num_stencil_x)
u_array = np.where(((x_array >= 30) |(x_array < 10)), 0.0, 1.0)
u_lower_boundary = 0.0
u_upper_boundary = 0.0
Time_step = 200
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
total_movement = C * Delta_t * (Time_step+1)
exact_u_array = np.where(((x_array >= 30 + total_movement) |(x_array < 10 + total_movement)), 0.0, 1.0)
plt.plot(x_array, u_array, label="Initial condition")
plt.plot(x_array, exact_u_array, label="Answer")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

1. FTCS(Forward in Time and Central Difference in Space) $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) $

Die Lösung vibriert. Ich bekomme überhaupt keine Antwort

u_ftcs = u_array.copy()
#Zeitschritt einstellen
for n in range(Time_step):
    u_old = u_ftcs.copy()
    u_ftcs[0] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary)
    u_ftcs[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1])
    for j in range(1,Num_stencil_x-1, 1):
        u_ftcs[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_ftcs, label="FTCS")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(FTCS)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

2. Differenzierung gegen den Wind

Die numerische Diffusion ist groß und die Lösung wird glatt.

u_upwind = u_array.copy()
#Zeitschritt einstellen
for n in range(Time_step):
    u_old = u_upwind.copy()
    u_upwind[0:1] = u_old[0] - CFL * (u_old[1] - u_lower_boundary)
    for j in range(1,Num_stencil_x):
        u_upwind[j:j+1] = u_old[j] - CFL * (u_old[j] - u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_upwind, label="Upwind differencing")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Upwind differencing)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

3. Lax-Wendroff scheme $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) + \frac{c^2}{2} \left( \frac{\Delta t}{\Delta x} \right)^2 \left(u_{j+1}^n - 2 u_j^n+ u_{j-1}^n \right) $

Es gibt numerische Vibration und stumpfe Lösung.

u_lw= u_array.copy()
#Zeitschritt einstellen
for n in range(Time_step):
    u_old = u_lw.copy()
    u_lw[0:1] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary) \
                    + CFL**2 / 2 * (u_old[1] - 2 * u_old[0] + u_lower_boundary)
    u_lw[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1]) \
                    + CFL**2 / 2 * (u_upper_boundary - 2 * u_old[-1] + u_old[-2])
    for j in range(1,Num_stencil_x-1, 1):
        u_lw[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1]) \
                    + CFL**2 / 2 * (u_old[j+1] - 2 * u_old[j] + u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_lw, label="Lax-Wendroff scheme")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Lax-Wendroff scheme)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

4. CIP(Constrained Interpolation Profile scheme) method

Die Grundidee der CIP-Methode besteht darin, dass beim ersten Diskretisieren einer Funktion f die physikalischen Größen eines Punktes $ x_j $ und eines Punktes $ x_ {j-1} $ daneben reibungslos miteinander verbunden werden. Auf diese Weise können die Differenzwerte der jeweiligen physikalischen Größen $ f_j $ und $ f_ {j-1} $ für die Zeitentwicklung verwendet werden, ohne dass die Neigung zwischen den Gittern divergiert.

Mit anderen Worten, unter der Bedingung, dass der ** Funktionswert und der Differenzwert an den Gitterpunkten ** stetig sind, kann die Beziehung zwischen den beiden Gitterpunkten durch die kubische Komplementärfunktion $ F (x) $ ausgedrückt werden.

Im Bereich von $ x_ {j-1} <= x <x_j $ $ F(x) = a X^3 + b X^2 +c X + d \quad ,where \quad X = x - x_j $ Wenn $ F (x) $ durch Vervollständigung dritter Ordnung wie in definiert ist, ** unter der Bedingung, dass der Funktionswert und der Differenzwert auf dem Gitterpunkt stetig sind **,

Ist zufrieden, wenn Sie dies in die Formel von $ F (x) $ einsetzen und $ a, b, c, d $ finden,

a = \frac{\partial f_j^n / \partial x + \partial f_{j+1}^n / \partial x }{\Delta x^2} + \frac{2 (f_j^n - f_{j+1}^n)}{\Delta x^3} \\
b =  \frac{3 (f_{j+1}^n - f_j^n)}{\Delta x^2} - \frac{2( \partial f_j^n / \partial x + \partial f_{j+1}^n / \partial x )}{\Delta x} \\
c = \frac{\partial f_j^n}{\partial x}\\
d = f_j^n

Wenn die Geschwindigkeit c konstant ist, erfüllt $ \ frac {\ partielles u} {\ partielles t} + c \ frac {\ partielles u} {\ partielles x} = 0 $ auch den Differenzwert, also die Zeitschrittbreite $ \ Delta t (n + 1 Schritt) Der Wert und der Differenzwert nach $ Sekunden sind

u_j^{n+1} = F(x_j - c \Delta t)= a \xi^3 + b \xi^2 + c\xi + d \\
\frac{\partial u_j^{n+1}}{\partial x} = \frac{d F(x_j - c \Delta t)}{d x} = 3 a \xi^2 + 2 b \xi + c \\
, where \quad \xi = - c \Delta t

Die numerische Diffusion und Vibration sind überwiegend kleiner als die anderen drei Methoden. ~~ Die Übertragung war jedoch nicht so gut wie Referenzen. ~~ (2020/02/20 Nachdem ich die Korrekturen vorgenommen habe, die ich in den Kommentaren erhalten habe, hat es sich erheblich verbessert.)

u_cip= u_array.copy()
partial_u_cip = ((np.append(u_cip[1:], u_upper_boundary) + u_cip)/2 - (np.append(u_lower_boundary, u_cip[:-1]) + u_cip)/2)/ Delta_x
#Zeit entwickeln
for n in range(Time_step):
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = 0
    partial_u_cip[0] = 0
    for j in range(1, Num_stencil_x):
        a = (partial_u_old[j] + partial_u_old[j-1]) / Delta_x**2 - 2.0 * (u_old[j] - u_old[j-1]) / Delta_x**3
        b = 3 * (u_old[j-1] - u_cip[j]) / Delta_x**2 + (2.0*partial_u_old[j] + partial_u_old[j-1]) / Delta_x
        c = partial_u_old[j]
        d = u_old[j]
        xi = - C * Delta_t  # C > 0
        u_cip[j:j+1] = a * xi**3 + b * xi**2 + c * xi + d
        partial_u_cip[j:j+1] = 3 * a * xi**2 + 2 * b * xi + c
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_cip, label="CIP method")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(CIP)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

Zusammenfassung

Die Grundlagen der Translokationsgleichungen, aus denen die Fluidgleichung besteht, werden mit einer kurzen Implementierung zusammengefasst. Die folgenden vier wurden implementiert.

• FTCS(Forward in Time and Central Difference in Space)
• Gegenwinddifferenzierung
• Lax-Wendroff scheme
• CIP(Constrained Interpolation Profile scheme) method

Nächstes Mal werde ich über [Diffusionsgleichung] sprechen (https://qiita.com/KQTS/items/97daa509991bb9777a4a).

Verweise

1. http://zellij.hatenablog.com/entry/20140805/p1
2. https://www.cradle.co.jp/media/column/a233
3. http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%B0%DC%CE%AE%CB%A1
4. http://www.astro.phys.s.chiba-u.ac.jp/netlab/mhd_summer2001/cip-intro.pdf
5. http://www.astro.phys.s.chiba-u.ac.jp/netlab/summer-school_2004/TEXTBOOK/text3-1.pdf