[Python] Fluidsimulation: Von linear zu nichtlinear

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. Ich bin der Meinung, dass die numerische Fluiddynamik eine ziemlich schwierige Studie ist, daher möchte ich sie so schreiben, dass sie für Anfänger leicht verständlich ist. 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

Serie

Grober Inhalt dieses Artikels

Als ersten Schritt zur Erstellung der für die Wassersimulation erforderlichen Grundgleichung der Flüssigkeit wird der Umgang mit nichtlinearen Gleichungen erläutert. Schließlich werden wir eine numerische Simulation der Bergers-Gleichung durchführen, die eine nichtlineare Gleichung ist.

Verwendete Bibliothek: Numpy, Scipy, Matplotlib

Ich werde diese zweidimensionale Version mit so einem Kerl machen.

burgers.gif

Inhaltsverzeichnis

Kapitel Titel Bemerkungen
1. Linear und nichtlinear
1-2. Merkmale linearer und nichtlinearer Gleichungen
1-3. Unterscheidung zwischen linear und nichtlinear Bean kenntnisreich
2. Implementierung der linearen Transferdiffusionsgleichung
2-1. Lineare Transferdiffusionsgleichung
2-2. Implementierung
3. Implementierung der Bergers-Gleichung
3-1. Was ist die Burger-Gleichung?? Entsprechung zur nichtlinearen Übertragungsdiffusionsgleichung
3-2. Implementierung
3-3. Ergänzung zum diskreten Ausdruck
4. Auf 2 Dimensionen erweitern
4-1. Zweidimensionale Übertragungsgleichung
4-2. Zweidimensionale Diffusionsgleichung
4-3. Zweidimensionale Burger-Gleichung

1. Linear und nichtlinear

1-1. Eigenschaften linearer und nichtlinearer Gleichungen

Das Folgende ist eine kurze Zusammenfassung von linearen und nichtlinearen Gleichungen.

Gleichung Charakteristisch
linear 1.Es gilt eine proportionale Beziehung
2. Kann eine Lösung finden
nicht linear 1.Proportionalbeziehung gilt nicht
2. Grundsätzlich kann ich keine Lösung finden(Es können nur wenige Gleichungen wie die kdV-Gleichung gelöst werden)

Flüssigkeitsgleichungen, die in Zukunft behandelt werden sollen, werden als nichtlineare Gleichungen klassifiziert, und im Grunde kann keine Lösung erhalten werden. Wenn Sie die Lösung einer nichtlinearen Gleichung finden, besteht die Grundrichtlinie darin, eine numerische Berechnung unter Verwendung der in den vorherigen Artikeln eingeführten diskreten Gleichungen durchzuführen und die Lösung ungefähr zu finden. ..

1-2. Wie man zwischen linear und nichtlinear unterscheidet

Neben dem Thema möchte ich erklären, wie man zwischen linearen und nichtlinearen Gleichungen unterscheidet. Um zwischen linearen und nichtlinearen Gleichungen zu unterscheiden, bestimmen Sie, ob die Lösungsüberlagerung gilt.

Ich möchte dies anhand der in Vorheriger Artikel eingeführten Translokationsgleichung als Beispiel erläutern.

Die Advance-Gleichung lautet $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $ Es wurde vertreten durch. Die Übertragungsrate c ist jedoch eine Konstante (was bedeutet, dass sich die Substanz mit einer konstanten Rate von c bewegt).

Wenn diese Gleichung zwei Lösungen $ u_1 $ und $ u_2 $ enthält

\frac{\partial u_1}{\partial t} + c \frac{\partial u_1}{\partial x} = 0
\frac{\partial u_2}{\partial t} + c \frac{\partial u_2}{\partial x} = 0

Ist festgelegt. Daher gilt für $ u_1 + u_2 , was die Summe der beiden Lösungen ist, $ \frac{\partial (u_1+u_2)}{\partial t} + c \frac{\partial (u_1+u_2)}{\partial x} = 0 $$ Ist etabliert, und wir können sehen, dass diese Gleichung linear ist.

Ändern wir nun die konstante Übertragungsgeschwindigkeit c in die Variable u.

\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0
\frac{\partial (u_1+u_2)}{\partial t} + (u_1+u_2) \frac{\partial (u_1+u_2)}{\partial x} \neq \frac{\partial (u_1+u_2)}{\partial t} + u_1 \frac{\partial u_1}{\partial x}+ u_2 \frac{\partial u_2}{\partial x} = 0

Dann ist in dieser Gleichung $ u_1 + u_2 $, die Summe der beiden Lösungen, keine Lösung, daher wissen wir, dass es sich um eine nichtlineare Gleichung handelt.

In diesem Artikel möchte ich mich mit solchen nichtlinearen Gleichungen befassen.

2. Implementierung der linearen Transferdiffusionsgleichung

2-1 Lineare Übertragungsdiffusionsgleichung

Zunächst möchte ich mich mit der linearen Transferdiffusionsgleichung befassen, die bisher auch als Entwicklungsbedeutung der Reihe dient. Diese Gleichung ist wie eine Kombination aus einer Übertragungsgleichung (eine Gleichung, die die Bewegung einer Substanz bedeutet) und einer Diffusionsgleichung (eine Gleichung, die die Gleichmäßigkeit physikalischer Größen bedeutet), wie in der folgenden Gleichung (c und) gezeigt. $ \ Nu $ ist eine Konstante). Die nichtlineare Version davon ist übrigens die Burgers-Gleichung, die im nächsten Kapitel erscheinen wird.

2-2. Implementierung

Vorerst werden wir es mit der CIP-Methode und der instationären iterativen Methode implementieren. Einfach ausgedrückt ist die CIP-Methode eine Methode zur Vorhersage zukünftiger Werte aus dem aktuellen Wert und seinem Differenzwert. Einzelheiten finden Sie im Artikel und in diesem pdf Siehe jp / netlab / mhd_summer2001 / cip-intro.pdf). Die instationäre iterative Methode ist eine Methode, um die Lösung eindimensionaler simultaner Gleichungen zu finden, während der Koeffizient geändert wird. Wenn Sie mehr wissen möchten, habe ich [Artikel über Diffusionsgleichungen] geschrieben (https://qiita.com/KQTS). Siehe / items / 97daa509991bb9777a4a) und Artikel über die instationäre Iterationsbibliothek von Python.

Es ist nicht so schwer zu tun. Bei der CIP-Methode werden die Translokationsphase und die Nicht-Translokationsphase getrennt betrachtet. Kurz gesagt, Sie müssen nur u aus der Translokationsgleichung berechnen und das Berechnungsergebnis u verwenden, um die Diffusionsgleichung zu lösen. In Bezug auf das Bild sind die Geschwindigkeit des Reisens durch Translokation und die Geschwindigkeit der Vereinheitlichung durch Diffusion unterschiedlich? (Frequenz?) Ich denke, es ist, als würde man getrennt denken. Die Translokationsdiffusionsgleichung kann wie folgt ausgedrückt werden, indem sie in eine Translokationsphase und eine Nicht-Translokationsphase unterteilt wird. Bei Verwendung der CIP-Methode sind auch Informationen zum Differenzwert $ \ frac {\ partielles u} {\ partielles x} $ erforderlich, sodass diese gleichzeitig mit der Berechnung von u durchgeführt werden.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \\
    \frac{\partial (\partial u / \partial x)}{\partial t} + c \frac{\partial (\partial u / \partial x)}{\partial x} = 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t}  = \nu \frac{\partial^2 u}{\partial x^2} \\
    \frac{\partial (\partial u / \partial x)}{\partial t} = \frac{\partial \left(\nu \frac{\partial^2 u}{\partial x^2}\right)}{\partial x}
 \end{array}
\right.

linear_adv_diff.gif

import sys
import os
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.sparse
import scipy.sparse.linalg as spla
from scipy.sparse import csr_matrix
from matplotlib.animation import ArtistAnimation
from mpl_toolkits.mplot3d import Axes3D

# Make stencils
# 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 = 300
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Nu = 0.5
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
diffusion_number = Nu * Delta_t / Delta_x**2
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.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)
print("Δx:", Delta_x, "Δt:", Delta_t, "CFL:", CFL, "ν:", Nu, "d:", diffusion_number)

def solve_advection_cip(u_cip, partial_u_cip, lower_boundary, velocity=C):
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = lower_boundary
    partial_u_cip[0] = 0  #Der Grenzwert wurde auf 0 gesetzt
    a = (partial_u_old[1:] + partial_u_old[:-1]) / Delta_x**2 - 2.0 * np.diff(u_old, n=1) / Delta_x**3
    b = 3 * (u_old[:-1] - u_cip[1:]) / Delta_x**2 + (2.0*partial_u_old[1:] + partial_u_old[:-1]) / Delta_x
    c = partial_u_old[1:]
    d = u_old[1:]
    xi = - velocity * Delta_t  # C > 0
    u_cip[1:] = a * xi**3 + b * xi**2 + c * xi + d
    partial_u_cip[1:] = 3 * a * xi**2 + 2 * b * xi + c
    return u_cip, partial_u_cip

def solve_diffusion(u_array, upper_boundary, lower_boundary, diffusion_number):
    a_matrix = np.identity(len(u_array)) * 2 *(1/diffusion_number+1) \
                - np.eye(len(u_array), k=1) \
                - np.eye(len(u_array), k=-1)
    temp_u_array = np.append(np.append(
                        lower_boundary, 
                        u_array), upper_boundary)
    b_array = 2 * (1/diffusion_number - 1) * u_array + temp_u_array[2:] + temp_u_array[:-2]
    b_array[0] += lower_boundary
    b_array[-1] += upper_boundary
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    u_array = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_array

fig = plt.figure()
images = []
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_cip_star, partial_u_cip_star = solve_advection_cip(u_cip, partial_u_cip, u_lower_boundary)
    u_cip = solve_diffusion(u_cip_star, u_upper_boundary, u_lower_boundary, diffusion_number)
    partial_u_cip = solve_diffusion(partial_u_cip_star, 0, 0, diffusion_number)  #Der Grenzwert wurde auf 0 gesetzt
    if n % 10 == 0:
        line, = plt.plot(x_array, u_cip, "r")
        images.append([line])
# plt.plot(x_array, u_array, label="Initial condition")
# plt.plot(x_array, u_cip, label="CIP method")
animatoin_gif = ArtistAnimation(fig, images)
animatoin_gif.save('linear_adv_diff.gif', writer="pillow")
plt.show()

3. Implementierung der Bergers-Gleichung

3-1. Wie lautet die Burgers-Gleichung?

Es ist eine der nichtlinearen partiellen Differentialgleichungen und eine nichtlineare Version der im vorherigen Kapitel erläuterten Translokationsdiffusionsgleichung. Wenn Sie Annahmen über die Fluidgleichung treffen und diese erheblich vereinfachen, wird sie zur Burgers-Gleichung. Diese Gleichung ist eine nichtlineare Gleichung, es ist jedoch bekannt, dass eine genaue Lösung unter Verwendung einer Call-Hop-Transformation erhalten werden kann. Daher wird es häufig verwendet, um die Gültigkeit numerischer Berechnungsmethoden zu untersuchen.

3-2. Implementierung

Implementiert die Burgers-Gleichung. Dies wird auch durch die CIP-Methode implementiert. Es ist wie folgt, wenn es in eine Translokationsphase und eine Nicht-Translokationsphase unterteilt wird.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 \\
    \frac{\partial (\partial u / \partial x)}{\partial t} + u \frac{\partial (\partial u / \partial x)}{\partial x} = 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t}  = \nu \frac{\partial^2 u}{\partial x^2} \\
    \frac{\partial (\partial u / \partial x)}{\partial t} = \frac{\partial \left(\nu \frac{\partial^2 u}{\partial x^2}\right)}{\partial x} - \left(\frac{\partial u}{\partial x} \right)^2
 \end{array}
\right.

Es gibt nicht viele Änderungen gegenüber der linearen Übertragungsdiffusionsgleichung. Wo man sich ändert

Ich denke darüber nach. Wir haben einige Änderungen an der Verbreitung vorgenommen. Weitere Informationen finden Sie im nächsten Abschnitt.

Wenn Sie die Berechnungsergebnisse kurz betrachten, können Sie sehen, dass die Geschwindigkeit umso höher ist, je höher die Höhe ist. Die obere Seite der Welle holt die untere Seite auf dem Weg ein und bildet eine steile Diskontinuitätsfläche ähnlich einer Stoßwelle. Danach können Sie sehen, dass die Lösung aufgrund des Diffusionseffekts allmählich stumpf wird.

burgers.gif

def solve_advection_cip(u_cip, partial_u_cip, lower_boundary, dt=Delta_t, velocity=C):
    """
Lösen der Übertragungsgleichung nach der CIP-Methode
    """
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = lower_boundary
    partial_u_cip[0] = (u_cip[0] - lower_boundary) / Delta_x
    a = (partial_u_old[1:] + partial_u_old[:-1]) / Delta_x**2 - 2.0 * np.diff(u_old, n=1) / Delta_x**3
    b = 3 * (u_old[:-1] - u_cip[1:]) / Delta_x**2 + (2.0*partial_u_old[1:] + partial_u_old[:-1]) / Delta_x
    c = partial_u_old[1:]
    d = u_old[1:]
    xi = - velocity * dt  # C > 0
    u_cip[1:] = a * xi**3 + b * xi**2 + c * xi + d
    partial_u_cip[1:] = 3 * a * xi**2 + 2 * b * xi + c
    return u_cip, partial_u_cip

def solve_diffusion(u_array, lower_boundary, upper_boundary, diffusion_number, 
                    update_partial=False, u_array_integral=None, prop_lambda=1/2):
    """
Lösen Sie die Diffusionsgleichung.
    
    Parameters
    ----------
    u_array : array-like
n Vektor
    upper_boundary : numpy.float64
        u_array[n]Neben an
    lower_boundary : numpy.float64
        u_array[0]Neben an
    diffusion_number : numpy.float64
Nummer verbreiten. Positive Zahl.
    update_partial : Bool, default False
Beim Aktualisieren der Differentialformel auf True setzen.
    u_array_integral : array-like, default None
Wird beim Aktualisieren der Differentialformel verwendet.[lower_boundary, u_array, upper_boudary]N in Form von+Vektor von 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Semi-implizite Methode. Kurbel-Nicolson scheme
        1: Fully implicit method.Zeitrückzugsunterschied.

    Returns
    -------
    u_array : array-like
n-reihige Matrix
    """
    a_matrix = np.identity(len(u_array)) * (1/diffusion_number+2 * prop_lambda) \
                - prop_lambda * np.eye(len(u_array), k=1) \
                - prop_lambda * np.eye(len(u_array), k=-1)
    temp_u_array = np.append(np.append(
                             lower_boundary, 
                             u_array), upper_boundary)
    b_array = (1/diffusion_number - 2 * (1-prop_lambda)) * u_array + (1-prop_lambda)*temp_u_array[2:] + (1-prop_lambda)*temp_u_array[:-2]
    b_array[0] += prop_lambda * lower_boundary
    b_array[-1] += prop_lambda * upper_boundary
    if update_partial:
        #Tun Sie dies, wenn Sie das Differential aktualisieren
        b_array -= ((u_array_integral[2:] - u_array_integral[:-2]) / 2 / Delta_x)**2
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    u_array = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_array

fig = plt.figure(1)
images = []
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
#Stabilitätsbedingungen einstellen
cfl_condition = 1
diffusion_number_restriction = 1/2
#Zeit entwickeln
for n in range(Time_step):
    delta_t = Delta_t
    cfl_max = np.max(np.abs(u_cip)) * delta_t / Delta_x
    diffusion_max = Nu * delta_t / Delta_x**2
    if cfl_max > cfl_condition:
        #Beurteilung der CFL-Bedingungen
        # cfl_Wenn es größer als die Bedingung ist, reduzieren Sie die Zeitschrittbreite dt, damit die CFL-Bedingung erfüllt ist.
        delta_t = cfl_condition * Delta_x / np.max(np.abs(u_cip))
    if diffusion_number > diffusion_number_restriction:
        #Beurteilung der Diffusionszahl
        # diffusion_number_Wenn es größer als die Beschränkung ist, verringern Sie die Zeitschrittbreite dt, so dass die Bedingung der Diffusionszahl erfüllt ist.
        delta_t = diffusion_max * Delta_x ** 2 / Nu
    #Lösen Sie die Translokationsgleichung
    u_cip_star, partial_u_cip_star = solve_advection_cip(u_cip, partial_u_cip, u_lower_boundary, delta_t, u_cip[1:])
    #Lösen Sie die Diffusionsgleichung
    u_cip = solve_diffusion(u_cip_star, u_lower_boundary, u_upper_boundary, diffusion_number)
    u_cip_with_boundary = np.append(np.append(
                                     u_lower_boundary, 
                                     u_cip_star), 
                                     u_upper_boundary)
    partial_u_cip = solve_diffusion(partial_u_cip_star, 0, 0, diffusion_number, 
                                    update_partial=True, 
                                    u_array_integral=u_cip_with_boundary)  #Der Grenzwert wurde auf 0 gesetzt
    if n % 10 == 0:
        line, = plt.plot(x_array, u_cip, "r")
        images.append([line])
# plt.plot(x_array, u_array, label="Initial condition")
# plt.plot(x_array, u_cip, label="Burgers")
animatoin_gif = ArtistAnimation(fig, images)
animatoin_gif.save('burgers.gif', writer="pillow")

3-3 Diskretionsformel der Diffusionsphase

Diffusionsgleichung

  \frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2}

Wird durch einen diskreten Ausdruck ausgedrückt

\frac{u_j^{n+1} - u_j^n}{\Delta t} = \nu \left((1 - \lambda) \frac{u_{j+1}^n - 2 u_j^n + u_{j-1}^n }{\Delta x^2} + \lambda \frac{u_{j+1}^{n+1} - 2 u_j^{n+1} + u_{j-1}^{n+1} }{\Delta x^2}\right)

Es wird sein. In Anbetracht der Zukunft wird der diskrete Ausdruck im Gegensatz zu Artikel zur Diffusionsgleichung mit $ \ lambda $ ausgedrückt. Um den Wert dieses $ \ lambda $ zusammenzufassen

\lambda Diskriminierungsmethode Charakteristisch
0 Mittenunterschied Explizite Methode. Berechnen Sie die Zeitdifferenz nur mit dem aktuellen Zeitwert.
1/2 Crank Nicholsons implizite Methode Semi-implizite Methode. Berechnen Sie die Zeitdifferenz mit der Hälfte des aktuellen Werts und der Hälfte des Werts beim nächsten Mal.
1 Zeitrückzugsunterschied Vollständig implizite Methode. Berechnen Sie die Zeitdifferenz nur mit dem nächsten Zeitwert.

Es wird sein. Unter dem Strich können Sie mit $ \ lambda $ drei verschiedene Verteilungstechniken implementieren.

Wenn Sie dies transformieren und den Wert der Zeit n + 1 auf die linke Seite bringen $ -\lambda u_{j+1}^{n+1} +\left(\frac{1}{d} + 2 \lambda \right) u_j^{n+1} - \lambda u_{j-1}^{n+1} = (1-\lambda) u_{j+1}^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_j^{n} + (1 - \lambda) u_{j-1}^{n} \\ d = \nu \frac{\Delta t}{\Delta x^2} $

Unter der Annahme, dass die Gitterpunkte im Bereich von 1 bis M liegen, wird der Grenzwert durch $ u_0, u_ {M + 1} $ dargestellt und in eine Matrixanzeige konvertiert.

  \left(
    \begin{array}{cccc}
      \frac{1}{d} + 2 \lambda & -\lambda                      &  0                      & \ldots   & \ldots                   & 0 \\
      -\lambda                & \frac{1}{d} + 2 \lambda       & -\lambda                & 0        & \ldots                   & 0 \\
      0                       &-\lambda                       & \frac{1}{d} + 2 \lambda & -\lambda & 0                         & \ldots  \\
      \vdots                  & \ddots                        & \ddots                  & \ddots   & \ddots                   & \ddots\\
      0                       & 0                             & \ddots                  & -\lambda & \frac{1}{d} + 2 \lambda & -\lambda \\
      0                       & 0                             & \ldots                  & 0        & -\lambda                 & \frac{1}{d} + 2 \lambda
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      u_1^{n+1}  \\
      u_2^{n+1}  \\
      u_3^{n+1}    \\
      \vdots \\
      u_{M-1}^{n+1} \\
      u_M^{n+1} 
    \end{array}
  \right)

   = 
   \left(
    \begin{array}{c}
      (1 - \lambda) u_2^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_1^{n} + \left((1-\lambda) u_0^n + \lambda u_0^{n+1}\right) \\
      (1-\lambda) u_3^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_2^{n} + (1-\lambda) u_1^n  \\
      (1-\lambda) u_4^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_3^{n} + (1-\lambda)u_2^n  \\
      \vdots \\
      (1-\lambda)u_M^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_{M-1}^{n} + (1-\lambda)u_{M-2}^n  \\
      \left((1-\lambda)u_{M+1}^n + \lambda u_{M+1}^{n+1}\right) + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_{M}^{n} + (1-\lambda)u_{M-1}^n
    \end{array}
  \right)

Wenn Sie $ \ frac {\ partielles u} {\ partielles x} $ aktualisieren möchten, ändern Sie u in dieser Matrix in $ \ frac {\ partielles u} {\ partielles x} $ und verschieben Sie es nach rechts.

\left(
    \begin{array}{c}
      - \left(\frac{u_2^{n} - u_0^{n}}{2 \Delta x} \right)^2  \\
      - \left(\frac{u_3^{n} - u_1^{n}}{2 \Delta x} \right)^2  \\
      - \left(\frac{u_4^{n} - u_2^{n}}{2 \Delta x} \right)^2    \\
      \vdots \\
      - \left(\frac{u_M^{n} - u_{M-1}^{n}}{2 \Delta x} \right)^2 \\
      - \left(\frac{u_{M+1}^{n} - u_M^{n}}{2 \Delta x} \right)^2 
    \end{array}
  \right)

Hinzufügen.

Dies ist im vorherigen Abschnitt implementiert. Die BiCG STAB-Methode wurde als iterative Lösungsmethode verwendet. Der Grund ist, dass der Name cool ist.

4. Auf 2 Dimensionen erweitern

Ich möchte ein cooles Bild zeichnen, deshalb werde ich es in diesem Kapitel zweidimensional machen. Es ist nur eine Formelerweiterung.

4-1 Zweidimensionale Übertragungsgleichung

Wie im Beispiel werden wir es mit der CIP-Methode implementieren. Ausführliche Informationen zur Erweiterung der Formel finden Sie in der Referenz unter 2D-CIP-Methode. Da es 10 Unbekannte gibt, müssen wir nur 10 Gleichungen erstellen. Sie müssen zum Zeitpunkt der Differenzierung nicht an die (i + 1, j + 1) -Koordinaten denken. In den Referenzen habe ich es gut gemacht, um den Rechenaufwand zu reduzieren, aber es ist schwer zu verstehen, deshalb werde ich es hier ehrlich machen.

Genau wie in einer Dimension

F(X,Y) = a_3 X^3 + a_2 X^2 + a_1 X + b_3 Y^3 + b_2 Y^2 + b_1 Y + c_3 X^2 Y + c_2 X Y^2 + c_1 XY + d \quad ,where \quad X = x - x_j

Definieren Sie $ F (X, Y) $ mit kubischer Vervollständigung wie in. Aus der Bedingung, dass ** Funktionswert und Differenzwert an Gitterpunkten stetig sind **

F(0, 0) = f_{i,j}, \quad F(\Delta x, 0) = f_{i+1, j}, \quad F(0, \Delta y) = f_{i, j+1}, \quad F(\Delta x, \Delta y) = f_{i+1, j+1}\\ \frac{\partial F(0,0)}{\partial x} = \frac{\partial f_{i,j}}{\partial x}, \quad \frac{\partial F(\Delta x, 0)}{\partial x} = \frac{\partial f_{i+1, j}}{\partial x},\quad \frac{\partial F(0, \Delta y)}{\partial x} = \frac{\partial f_{i, j+1}}{\partial x} \\ \frac{\partial F(0,0)}{\partial y} = \frac{\partial f_{i,j}}{\partial y}, \quad \frac{\partial F(\Delta x, 0)}{\partial y} = \frac{\partial f_{i+1, j}}{\partial y},\quad \frac{\partial F(0, \Delta y)}{\partial y} = \frac{\partial f_{i, j+1}}{\partial y} \\ , where \quad \Delta x = x_{i+1} - x_i, \quad \Delta y = y_{j+1} - y_{j}

Die 10 Formeln von Der Wert am Gitterpunkt (i, j) der Funktion f und der Differenzwert werden jedoch ausgedrückt als $ f_ {i, j}, \ frac {\ partiell f_ {i, j}} {\ partiell x} $. Ich werde. Setzen Sie dies in die Formel von $ F (X, Y) $ ein, um den Koeffizienten zu finden. $ a_3 = \frac{\partial f_{i+1,j} / \partial x + \partial f_{i, j} / \partial x }{\Delta x^2} + \frac{2 (f_{i,j} - f_{i+1, j})}{\Delta x^3} \\ a_2 = \frac{3 (f_{i+1, j} - f_{i,j})}{\Delta x^2} - \frac{(2 \partial f_{i,j} / \partial x + \partial f_{i+1, j} / \partial x )}{\Delta x} \\ a_1 = \frac{\partial f_{i,j}}{\partial x}\\ b_3 = \frac{\partial f_{i,j+1} / \partial y + \partial f_{i, j} / \partial y }{\Delta y^2} + \frac{2 (f_{i,j} - f_{i, j+1})}{\Delta y^3} \\ b_2 = \frac{3 (f_{i, j+1} - f_{i,j})}{\Delta y^2} - \frac{(2 \partial f_{i,j} / \partial y + \partial f_{i, j+1} / \partial y )}{\Delta y} \\ b_1 = \frac{\partial f_{i,j}}{\partial y}\\ c_3 = \frac{\partial f_{i+1,j}/\partial y - \partial f_{i,j}/\partial y }{(\Delta x)^2} - \frac{c_1}{\Delta x}\\ c_2 = \frac{\partial f_{i,j+1}/\partial x - \partial f_{i,j}/\partial x }{(\Delta y)^2} - \frac{c_1}{\Delta y}\\ c_1 = - \frac{f_{i+1,j+1} - f_{i, j+1} - f_{i+1,j} + f_{i, j}}{\Delta x \Delta y} + \frac{\partial f_{i+1,j}/\partial y - \partial f_{i,j}/\partial y }{\Delta x} + \frac{\partial f_{i,j+1}/\partial x - \partial f_{i,j}/\partial x }{\Delta y}\\ d = f_{i,j} $

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

    u_{i,j}^{n+1} = F(x_i - c_1 \Delta t, y_j -c_2 \Delta t)= a_3 \xi^3 + a_2 \xi^2 + a_1 \xi + b_3 \zeta^3 + b_2 \zeta^2 + b_1 \zeta + c_3 \xi^2 \zeta + c_2 \xi \zeta^2 + c_1 \xi \zeta + d \\
    \frac{\partial u_{i,j}^{n+1}}{\partial x} = \frac{\partial F(x_i - c_1 \Delta t, y_j - c_2 \Delta t)}{\partial x} = 3 a_3 \xi^2 + 2 a_2 \xi + a_1 + 2 c_3 \xi \zeta c_2 \zeta^2 + c_1 \zeta \\
    \frac{\partial u_{i,j}^{n+1}}{\partial y} = \frac{\partial F(x_i - c_1 \Delta t, y_j - c_2 \Delta t)}{\partial y} = 3 b_3 \zeta^2 + 2 b_2 \zeta + b_1 + 2 c_2 \xi \zeta c_3 \xi^2 + c_1 \xi \\
    , where \quad \xi = - c_1 \Delta t, \quad \zeta = - c_2 \Delta t \quad (c_1, c_2 < 0)

Wenn die Geschwindigkeit nicht konstant ist $ \ frac {\ partielles u} {\ partielles t} + u \ frac {\ partielles u} {\ partielles x} + v \ frac {\ partielles u} {\ partielles y} = 0 $ Zusätzliche Begriffe, die nach der räumlichen Differenzierung erscheinen, werden als Nicht-Translokationsbegriffe berechnet.

2d-advection.gif

4-2 Zweidimensionale Diffusionsgleichung

Zweidimensionale Diffusionsgleichung

    \frac{\partial u}{\partial t} = \nu \left\{ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right\}

Wird durch einen diskreten Ausdruck ausgedrückt

  \frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} = \nu \left\{(1 - \lambda) \left( \frac{u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n }{\Delta x^2} + \frac{u_{i,j+1}^n - 2 u_{i,j}^n + u_{i,j-1}^n }{\Delta y^2} \right) + \lambda \left( \frac{u_{i+1,j}^{n+1} - 2 u_{i,j}^{n+1} + u_{i-1,j}^{n+1} }{\Delta x^2} + \frac{u_{i,j+1}^{n+1} - 2 u_{i,j}^{n+1} + u_{i,j-1}^{n+1} }{\Delta y^2} \right) \right\}

Sortieren Sie die Raster in einer Reihe, um die Implementierung in Ihren Code zu vereinfachen. Setzen Sie den Gitterpunkt des gesamten Systems auf $ M = NX \ times NY $ und setzen Sie die Position des Gitterpunkts (i, j) auf den m-ten (Start 0 gemäß Python-Notation), der vom Gitterpunkt (0,0) zählt. Ich werde. Dann ist der Gitterpunkt (i + 1, j) das m + 1 und der Gitterpunkt (i, j + 1) das m + NX. Wenden Sie dies auf die obige Formel an, wenn Sie den Wert der Zeit n + 1 auf die linke Seite bringen

- d_2 \lambda u_{m+NX}^{n+1} - d_1 \lambda u_{m+1}^{n+1}  +s u_{m}^{n+1} - d_1 \lambda u_{m-1}^{n+1} - d_2 \lambda u_{m-NX}^{n+1} = d_2 (1-\lambda) u_{m+NX}^{n} + d_1 (1-\lambda) u_{m+1}^{n} + t u_m^{n} + d_1 (1 - \lambda) u_{m-1}^{n} + d_2 (1 - \lambda) u_{m-NX}^{n} \\
,where \quad s = 1 + 2 \lambda (d_1 + d_2), \quad t = 1 - 2 (1-\lambda) (d_1 + d_2)\\
d_1 = \nu \frac{\Delta t}{\Delta x^2}, \quad d_2 = \nu \frac{\Delta t}{\Delta y^2}

Wenn Sie es also einschließlich des Grenzwerts implementieren, sind Sie fertig. Da die Matrixanzeige lang ist, werde ich sie hier überspringen.

2d-diffusion.gif

4-3 Zweidimensionale Burger-Gleichung

Machen wir es vorerst zweidimensional. Ändern Sie die Notation und u in f, um von der Geschwindigkeit zu unterscheiden.

\left\{
\begin{array}{l}
    \frac{\partial f}{\partial t} + u \frac{\partial f}{\partial x} + v \frac{\partial f}{\partial y} = 0 \\
    \frac{\partial (\partial f / \partial x)}{\partial t} + u \frac{\partial (\partial f / \partial x)}{\partial x} + v \frac{\partial (\partial f / \partial x)}{\partial x}= 0  \\
    \frac{\partial (\partial f / \partial y)}{\partial t} + u \frac{\partial (\partial f / \partial y)}{\partial x} + v \frac{\partial (\partial f / \partial y)}{\partial x}= 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial f}{\partial t}  = \nu \left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\right) \\
    \frac{\partial (\partial f / \partial x)}{\partial t} = \frac{\partial \left\{ \nu \left( \frac{\partial^2 f}{\partial x^2}　+ \frac{\partial^2 f}{\partial y^2}\right)\right\}}{\partial x} - \left(\frac{\partial u}{\partial x} \right)\left(\frac{\partial f}{\partial x} \right) - \left(\frac{\partial v}{\partial x} \right) \left(\frac{\partial f}{\partial y} \right)  \\
      \frac{\partial (\partial f / \partial y)}{\partial t} = \frac{\partial \left\{ \nu \left( \frac{\partial^2 f}{\partial x^2}　+ \frac{\partial^2 f}{\partial y^2}\right)\right\}}{\partial y} - \left(\frac{\partial u}{\partial y} \right)\left(\frac{\partial f}{\partial x} \right) - \left(\frac{\partial v}{\partial y} \right) \left(\frac{\partial f}{\partial y} \right) 
\end{array}
\right.

Implementieren Sie dies einfach. Wie unten gezeigt, weiß ich nicht, ob es übereinstimmt, aber Sie können ein solches Berechnungsergebnis erhalten.

2d-burgers_short.gif

# Make stencils
# Creat square wave
Num_stencil_x = 101
Num_stencil_y = 101
x_stencil = np.arange(Num_stencil_x)
y_stencil = np.arange(Num_stencil_y)
x_array, y_array = np.meshgrid(x_stencil, y_stencil)
u_array = np.where(((30 >x_array) & (x_array >= 10) & (30 > y_array) & (y_array >= 10)), 1.0, 0.0)
u_x_lower_boundary = np.zeros(Num_stencil_x)
u_y_lower_boundary = np.zeros(Num_stencil_y)
u_x_upper_boundary = np.zeros(Num_stencil_x)
u_y_upper_boundary = np.zeros(Num_stencil_y)
ux_x_lower_boundary = np.zeros(Num_stencil_x)
ux_y_lower_boundary = np.zeros(Num_stencil_y)
ux_x_upper_boundary = np.zeros(Num_stencil_x)
ux_y_upper_boundary = np.zeros(Num_stencil_y)
uy_x_lower_boundary = np.zeros(Num_stencil_x) 
uy_y_lower_boundary = np.zeros(Num_stencil_y)
uy_x_upper_boundary = np.zeros(Num_stencil_x)
uy_y_upper_boundary = np.zeros(Num_stencil_y)
coner_xlow_ylow = np.zeros(1)
coner_xup_ylow = np.zeros(1)
coner_xlow_yup = np.zeros(1)
coner_xup_yup = np.zeros(1)
Time_step = 300
Delta_x = max(x_stencil) / (Num_stencil_x-1)
Delta_y = max(y_stencil) / (Num_stencil_y-1)
C = 1
Nu = 0.5
Delta_t = 0.2
CFL_x = C * Delta_t / Delta_x
CFL_y = C * Delta_t / Delta_y
CFL = min(CFL_x, CFL_y)
diffusion_number_x = Nu * Delta_t / Delta_x**2
diffusion_number_y = Nu * Delta_t / Delta_y**2
total_movement = C * Delta_t * (Time_step+1)
print("Δx:", Delta_x, "Δt:", Delta_t, "CFL:", CFL_x, "ν:", Nu, "d:", diffusion_number_x)
fig = plt.figure(2)
axe = fig.add_subplot(111, projection='3d')
surface = axe.plot_surface(x_array, y_array, u_array, cmap="bwr", label="Initail condition")
fig.colorbar(surface)
axe.set_xlabel("x")
axe.set_ylabel("y")
axe.set_zlabel("u")
axe.set_xlim(0, 100)
axe.set_ylim(0, 100)
axe.set_zlim(0, 1)
axe.set_title("2d initial condition")
plt.show()

def solve_advection_cip_2d(u_cip, ux_cip, uy_cip, 
                           x_lower_boundary, x_upper_boundary, 
                           y_lower_boundary, y_upper_boundary,
                           ux_x_lower_boundary, ux_x_upper_boundary,
                           ux_y_lower_boundary, ux_y_upper_boundary,
                           uy_x_lower_boundary, uy_x_upper_boundary,
                           uy_y_lower_boundary, uy_y_upper_boundary,
                           coner_ll, coner_lu, coner_ul, coner_uu,
                           dt=Delta_t, 
                           velocity_x=C, velocity_y=0.0):
    """
    Solve the advection equation by using CIP method.
    
       -  x  +
     - 5333337
       1*****2
     y 1*****2
       1*****2
       1*****2
     + 6444448
     
    *: u_cip or ux_cip or uy_cip
    1: x_lower_boundary or ux_x_lower_boundary or uy_x_lower_boundary
    2: x_upper_boundary or ux_x_upper_boundary or uy_x_upper_boundary
    3: y_lower_boundary or ux_y_lower_boundary or uy_y_lowerboundary
    4: y_upper_boundary or ux_y_upper_boundary or uy_y_upper_boundary
    5: coner_ll
    6: coner_lu
    7: coner_ul
    8: coner_uu
    """
    if type(velocity_x) is not np.ndarray:
        velocity_x = np.ones(u_cip.shape) * velocity_x
    if type(velocity_y) is not np.ndarray:
        velocity_y = np.ones(u_cip.shape) * velocity_y
    # Memory the present values
    u_old = u_cip.copy()
    ux_old = ux_cip.copy()
    uy_old = uy_cip.copy()

    # Prepare for calculation
    xi = - np.sign(velocity_x) * velocity_x * dt
    zeta = - np.sign(velocity_y) * velocity_y * dt
    dx = - np.sign(velocity_x) * Delta_x
    dy = - np.sign(velocity_y) * Delta_y
    
    # Set the neighbourhood values for reducing the for loop
    u_with_xlower_bc = np.concatenate([x_lower_boundary[:, np.newaxis], u_old[:, :-1]], axis=1)
    u_with_xupper_bc = np.concatenate([u_old[:, 1:], x_upper_boundary[:, np.newaxis]], axis=1)
    u_with_ylower_bc = np.concatenate([y_lower_boundary[np.newaxis, :], u_old[:-1, :]], axis=0)
    u_with_yupper_bc = np.concatenate([u_old[1:, :], y_upper_boundary[np.newaxis, :]], axis=0)
    temp_u_with_all_bc = np.concatenate([x_lower_boundary[:, np.newaxis], u_old, 
                                         x_upper_boundary[:, np.newaxis]], axis=1)
    u_with_all_bc = np.concatenate([np.concatenate([coner_ll, y_lower_boundary,coner_ul])[np.newaxis, :],
                                    temp_u_with_all_bc,
                                    np.concatenate([coner_lu, y_upper_boundary, coner_uu])[np.newaxis, :]], axis=0)
    
    u_next_x = np.where(velocity_x > 0.0, u_with_xlower_bc, u_with_xupper_bc)
    u_next_y = np.where(velocity_y > 0.0, u_with_ylower_bc, u_with_yupper_bc)
    u_next_xy = np.where(((velocity_x > 0.0) & (velocity_y > 0.0)), u_with_all_bc[:-2, :-2], 
                    np.where(((velocity_x < 0.0) & (velocity_y > 0.0)), u_with_all_bc[:-2, 2:],
                        np.where(((velocity_x > 0.0) & (velocity_y < 0.0)), u_with_all_bc[2:, :-2], u_with_all_bc[2:, 2:])))
    ux_next_x = np.where(velocity_x > 0, 
                        np.concatenate([ux_x_lower_boundary[:, np.newaxis], ux_old[:, :-1]], axis=1), 
                        np.concatenate([ux_old[:, 1:], ux_x_upper_boundary[:, np.newaxis]], axis=1))
    ux_next_y = np.where(velocity_y > 0, 
                        np.concatenate([ux_y_lower_boundary[np.newaxis, :], ux_old[:-1, :]], axis=0), 
                        np.concatenate([ux_old[1:, :], ux_y_upper_boundary[np.newaxis, :]], axis=0))
    uy_next_x = np.where(velocity_x > 0, 
                        np.concatenate([uy_x_lower_boundary[:, np.newaxis], uy_old[:, :-1]], axis=1), 
                        np.concatenate([uy_old[:, 1:], uy_x_upper_boundary[:, np.newaxis]], axis=1))
    uy_next_y = np.where(velocity_y > 0, 
                        np.concatenate([uy_y_lower_boundary[np.newaxis, :], uy_old[:-1, :]], axis=0), 
                        np.concatenate([uy_old[1:, :], uy_y_upper_boundary[np.newaxis, :]], axis=0))
    u_next_x = np.where(velocity_x == 0.0, u_old, u_next_x)
    u_next_y = np.where(velocity_y == 0.0, u_old, u_next_y)
    u_next_xy = np.where(((velocity_x == 0.0) & (velocity_y != 0.0)), u_next_y, 
                    np.where(((velocity_x != 0.0) & (velocity_y == 0.0)), u_next_x,
                        np.where(((velocity_x == 0.0) & (velocity_y == 0.0)), u_old, u_next_xy)))
    ux_next_x = np.where(velocity_x == 0, ux_old, ux_next_x)
    ux_next_y = np.where(velocity_y == 0, ux_old, ux_next_y)
    uy_next_x = np.where(velocity_x == 0, uy_old, uy_next_x)
    uy_next_y = np.where(velocity_y == 0, uy_old, uy_next_y)
    dx = np.where(velocity_x == 0, 1.0, dx)
    dy = np.where(velocity_y == 0, 1.0, dy)
    
    # Calculate the cip coefficient
    c_1 = - (u_next_xy - u_next_x - u_next_y + u_old) / dx / dy \
            + (ux_next_y - ux_old) / dy + (uy_next_x - uy_old) / dx
    c_3 = (uy_next_x - uy_old) / dx ** 2 - c_1 / dx
    c_2 = (ux_next_y - ux_old) / dy ** 2 - c_1 / dy
    a_1 = ux_old
    a_2 = 3 * (u_next_x - u_old) / dx**2 - (ux_next_x + 2 * ux_old) / dx
    a_3 = (ux_next_x - ux_old) / 3 / dx**2 - 2 * a_2 / 3 / dx
    b_1 = uy_old
    b_2 = 3 * (u_next_y - u_old) / dy**2 - (uy_next_y + 2 * uy_old) / dy
    b_3 = (uy_next_y - uy_old) / 3 / dy**2 - 2 * b_2 / 3 / dy
    d = u_old
    
    # Calculate the values
    u_cip = a_3 * xi**3 + a_2 * xi**2 + a_1 * xi + \
            b_3 * zeta**3 + b_2 * zeta**2 + b_1 * zeta + \
            c_3 * xi**2 * zeta + c_2 * xi * zeta**2 + c_1 * xi * zeta + d
    ux_cip = 3 * a_3 * xi**2   + 2 * a_2 * xi + a_1 + 2 * c_3 * xi * zeta + c_2 * zeta**2 + c_1 * zeta
    uy_cip = 3 * b_3 * zeta**2 + 2 * b_2 * zeta + b_1 + 2 * c_2 * xi * zeta + c_2 * xi**2 + c_1 * xi
    return u_cip, ux_cip, uy_cip

def update_boundary_condition(u_cip, x_lower_boundary, x_upper_boundary, 
                             y_lower_boundary, y_upper_boundary):
    x_lower_boundary = min(0.0, np.min(u_cip[:, 0]))
    x_upper_boundary = max(0.0, np.max(u_cip[:, -1]))
    y_lower_boundary = min(0.0, np.min(u_cip[0, :]))
    y_upper_boundary = max(0.0, np.max(u_cip[-1, :]))

def solve_diffusion_2d(u_2d, 
                       x_lower_boundary, x_upper_boundary, 
                       y_lower_boundary, y_upper_boundary,
                       d_number_x, d_number_y,
                       v_x=None, v_y=None,
                       v_x_lower_bc=None, v_x_upper_bc=None,
                       v_y_lower_bc=None, v_y_upper_bc=None,
                       update_partial=0, u_array_integral=None, 
                       prop_lambda=1/2):
    """
Lösen Sie die Diffusionsgleichung. Es enthält auch Elemente von Nicht-Translokationsbegriffen.
    
    Parameters
    ----------
    u_2d : array-like
nx × ny Zeilenmatrix
    upper_boundary : numpy.float64
        u_array[n]Neben an
    lower_boundary : numpy.float64
        u_array[0]Neben an
    d_number_* : numpy.float64
Nummer verbreiten. Positive Zahl.
    update_partial : Bool, default False
Beim Aktualisieren der Differentialformel auf True setzen.
    u_array_integral : array-like, default None
Wird beim Aktualisieren der Differentialformel verwendet.[lower_boundary, u_array, upper_boudary]N in Form von+Vektor von 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Semi-implizite Methode. Kurbel-Nicolson scheme
        1: Fully implicit method.Zeitrückzugsunterschied.

    Returns
    -------
    u_2d : array-like
nx × ny Zeilenmatrix
    """
    #Holen Sie sich die Matrixgröße
    nx, ny = u_2d.shape
    matrix_size = nx * ny
    # Ax=Erstellen Sie A und b für b
    s_param = 1 + 2 * prop_lambda * (d_number_x + d_number_y)
    t_param = 1 - 2 * (1 - prop_lambda) * (d_number_x + d_number_y)
    east_matrix = np.eye(matrix_size, k=1)
    east_matrix[np.arange(nx-1, matrix_size-nx, nx), np.arange(nx, matrix_size-nx+1, nx)] *= 0
    west_matrix = np.eye(matrix_size, k=-1)
    west_matrix[np.arange(nx, matrix_size-nx+1, nx), np.arange(nx-1, matrix_size-nx, nx)] *= 0
    a_matrix = np.identity(matrix_size) * s_param \
                - prop_lambda * d_number_x * east_matrix \
                - prop_lambda * d_number_x * west_matrix \
                - prop_lambda * d_number_y * np.eye(matrix_size, k=nx) \
                - prop_lambda * d_number_y * np.eye(matrix_size, k=-nx)
    b_array = t_param * u_2d.flatten()
    temp_csr = d_number_x * (1 - prop_lambda) * (csr_matrix(east_matrix) + csr_matrix(west_matrix)) \
               + d_number_y * (1 - prop_lambda) * (csr_matrix(np.eye(matrix_size, k=nx)) + csr_matrix(np.eye(matrix_size, k=-nx)))
    b_array += temp_csr.dot(b_array)
    #Randbedingungen festlegen/ Setting of boundary condition
    b_array[:nx] += ((1 - prop_lambda) * y_lower_boundary + prop_lambda * y_lower_boundary) * d_number_y
    b_array[matrix_size-nx:] += ((1 - prop_lambda) * y_upper_boundary + prop_lambda * y_upper_boundary) * d_number_y
    b_array[np.arange(nx-1, matrix_size, nx)] += ((1 - prop_lambda) * x_upper_boundary + prop_lambda * x_upper_boundary) * d_number_x
    b_array[np.arange(0, matrix_size-nx+1, nx)] += ((1 - prop_lambda) * x_lower_boundary + prop_lambda * x_lower_boundary) * d_number_x
    if update_partial == 1:
        #Wird ausgeführt, wenn der Differentialausdruck in x-Richtung aktualisiert wird. Wird beim Lösen gewöhnlicher Diffusionsgleichungen nicht verwendet.
        if type(v_x) is not np.ndarray:
            v_x = np.ones(u_2d.shape) * v_x
        if type(v_y) is not np.ndarray:
            v_y = np.ones(u_2d.shape) * v_y
        v_x_with_bc = np.concatenate([v_x_lower_bc[:, np.newaxis], v_x, v_x_upper_bc[:, np.newaxis]], axis=1)
        v_y_with_bc = np.concatenate([v_y_lower_bc[:, np.newaxis], v_y, v_y_upper_bc[:, np.newaxis]], axis=1)
        u_with_xbc = np.concatenate([x_lower_boundary[:, np.newaxis], u_array_integral, x_upper_boundary[:, np.newaxis]], axis=1)
        u_with_ybc = np.concatenate([y_lower_boundary[np.newaxis, :], u_array_integral, y_upper_boundary[np.newaxis, :]], axis=0)
        temp_b_array = - ((v_x_with_bc[:, 2:] - v_x_with_bc[:, :-2]) / 2 / Delta_x) * \
                             ((u_with_xbc[:, 2:] - u_with_xbc[:, :-2]) / 2 / Delta_x)
        temp_b_array -= ((v_y_with_bc[:, 2:] - v_y_with_bc[:, :-2]) / 2 / Delta_x) * \
                               ((u_with_ybc[2:, :] - u_with_ybc[:-2, :]) / 2 / Delta_y)
        b_array += temp_b_array.flatten()
    elif update_partial == 2:
        #Wird ausgeführt, wenn der Differentialausdruck in y-Richtung aktualisiert wird. Wird beim Lösen gewöhnlicher Diffusionsgleichungen nicht verwendet.
        if type(v_x) is not np.ndarray:
            v_x = np.ones(u_2d.shape) * v_x
        if type(v_y) is not np.ndarray:
            v_y = np.ones(u_2d.shape) * v_y
        v_x_with_bc = np.concatenate([v_x_lower_bc[np.newaxis, :], v_x, v_x_upper_bc[np.newaxis, :]], axis=0)
        v_y_with_bc = np.concatenate([v_y_lower_bc[np.newaxis, :], v_y, v_y_upper_bc[np.newaxis, :]], axis=0)
        u_with_xbc = np.concatenate([x_lower_boundary[:, np.newaxis], u_array_integral, x_upper_boundary[:, np.newaxis]], axis=1)
        u_with_ybc = np.concatenate([y_lower_boundary[np.newaxis, :], u_array_integral, y_upper_boundary[np.newaxis, :]], axis=0)
        
        temp_b_array = - ((v_x_with_bc[2:, :] - v_x_with_bc[:-2, :]) / 2 / Delta_y) * \
                             ((u_with_xbc[:, 2:] - u_with_xbc[:, :-2]) / 2 / Delta_x) \
                       - ((v_y_with_bc[2:, :] - v_y_with_bc[:-2, :]) / 2 / Delta_y) * \
                               ((u_with_ybc[2:, :] - u_with_ybc[:-2, :]) / 2 / Delta_y)
        b_array += temp_b_array.flatten()
    a_csr = csr_matrix(a_matrix)
    u_2d = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_2d.reshape(nx, ny)

#2d Bugers Gleichung
images = []
fig_5 = plt.figure(5, figsize=(5,5))
axe = fig_5.add_subplot(111, projection='3d')
# surface = axe.plot_surface(x_array, y_array, u_array, cmap="bwr")
# fig.colorbar(surface)
axe.set_xlabel("x")
axe.set_ylabel("y")
axe.set_zlabel("u")
axe.set_xlim(0, 100)
axe.set_ylim(0, 100)
axe.set_zlim(0, 1)
axe.set_title("2d burgers")

u_cip = u_array.copy()
partial_u_cip_x \
    = (np.concatenate([u_cip[:, 1:], u_x_upper_boundary[:, np.newaxis]], axis=1) \
       - np.concatenate([u_x_lower_boundary[:, np.newaxis], u_cip[:, :-1]], axis=1)) / 2 / Delta_x
partial_u_cip_y \
    = (np.concatenate([u_cip[1:, :], u_y_upper_boundary[np.newaxis, :]], axis=0) \
       - np.concatenate([u_y_lower_boundary[np.newaxis, :], u_cip[:-1, :]], axis=0)) / 2 / Delta_y
u_boundary = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)

#Stabilitätsbedingungen einstellen
cfl_condition = 1
diffusion_number_restriction = 1/2
#Zeit entwickeln
for n in range(21):
    update_boundary_condition(u_cip, *u_boundary)
    delta_t = Delta_t
    cfl_max = np.max(np.abs(u_cip)) * delta_t / min(Delta_x, Delta_y)
    if cfl_max > cfl_condition:
        #Beurteilung der CFL-Bedingungen
        # cfl_Wenn es größer als die Bedingung ist, reduzieren Sie die Zeitschrittbreite dt, damit die CFL-Bedingung erfüllt ist.
        delta_t = CFL * min(Delta_x, Delta_y) / np.max(np.abs(u_cip))
    diffusion_max = Nu * delta_t / (Delta_x**2 + Delta_y**2)
    if diffusion_max > diffusion_number_restriction:
        #Beurteilung der Diffusionszahl
        # diffusion_number_Wenn es größer als die Beschränkung ist, verringern Sie die Zeitschrittbreite dt, so dass die Bedingung der Diffusionszahl erfüllt ist.
        delta_t = diffusion_max * (Delta_x ** 2 + Delta_y**2) / Nu
        diffusion_number_x *= delta_t / Delta_t
        diffusion_number_y *= delta_t / Delta_t
    u_props = (u_cip, partial_u_cip_x, partial_u_cip_y)
    u_boundary = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)
    ux_boundary = (ux_x_lower_boundary, ux_x_upper_boundary, ux_y_lower_boundary, ux_y_upper_boundary)
    uy_boundary = (uy_x_lower_boundary, uy_x_upper_boundary, uy_y_lower_boundary, uy_y_upper_boundary)
    u_coner = (coner_xlow_ylow, coner_xlow_yup, coner_xup_ylow, coner_xup_yup)
    diffusion_numbers = (diffusion_number_x, diffusion_number_y)
    u_cip_star, partial_u_cip_x_star, partial_u_cip_y_star \
        = solve_advection_cip_2d(*u_props, *u_boundary, *ux_boundary, *uy_boundary, *u_coner, dt=delta_t)
    velocities = (u_cip_star, 0.0)
    velocity_x_boundaries = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)
    velocity_y_boundaries = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary) #Vorerst
    u_cip = solve_diffusion_2d(u_cip_star, *u_boundary, *diffusion_numbers)
    partial_u_cip_x = solve_diffusion_2d(partial_u_cip_x_star, *ux_boundary, *diffusion_numbers,
                                         *velocities, *velocity_x_boundaries,
                                         update_partial=1, u_array_integral=u_cip_star)
    partial_u_cip_y = solve_diffusion_2d(partial_u_cip_y_star, *uy_boundary, *diffusion_numbers,
                                         *velocities, *velocity_y_boundaries,
                                         update_partial=2, u_array_integral=u_cip_star)
    
    if n % 1 == 0 and n > 0:
        img = axe.plot_wireframe(x_array, y_array, u_cip, cmap="bwr")
        images.append([img])
ani = animation.ArtistAnimation(fig_5, images, interval=100)
# ani.save('wf_anim_art.mp4',writer='ffmpeg',dpi=100)
ani.save('2d-burgers.gif', writer="pillow")
HTML(ani.to_html5_video())

Zusammenfassung

Das nächste Mal möchte ich mich mit Grundgleichungen für Flüssigkeiten befassen. Außerdem ist die Berechnung ziemlich umfangreich geworden, daher würde ich sie gerne beschleunigen, wenn ich Lust dazu habe.

Verweise

Recommended Posts

[Python] Fluidsimulation: Von linear zu nichtlinear
Änderungen von Python 3.0 zu Python 3.5
Änderungen von Python 2 zu Python 3.0
Post von Python nach Slack
Flirte von PHP nach Python
Anaconda aktualisiert von 4.2.0 auf 4.3.0 (python3.5 aktualisiert auf python3.6)
Wechseln Sie von Python2.7 zu Python3.6 (centos7)
Stellen Sie von Python aus eine Verbindung zu SQLite her
Rufen Sie Matlab von Python zur Optimierung auf
Post von Python auf Facebook Timeline
[Lambda] [Python] Von Lambda auf Twitter posten!
Stellen Sie von Python aus eine Verbindung zur utf8mb4-Datenbank her
Python (vom ersten Mal bis zur Ausführung)
Poste ein Bild von Python auf Tumblr
[Python] Flüssigkeitssimulation: Inkompressible Navier-Stokes-Gleichung
So greifen Sie über Python auf Wikipedia zu
Python, um von einer anderen Sprache zu wechseln
Hat sich nicht von Python 2 auf 3 geändert
Aktualisieren Sie Mac Python von 2 auf 3
Von Python bis zur Verwendung von MeCab (und CaboCha)
Einführung in die diskrete Ereignissimulation mit Python # 1
[Python] Flüssigkeitssimulation: Implementieren Sie die Übertragungsgleichung
So aktualisieren Sie Google Sheets von Python
Privates Python-Handbuch (von Zeit zu Zeit aktualisiert)
Ich möchte ein Glas aus Python verwenden
Konvertieren Sie von Katakana zu Vokal Kana [Python]
Push-Benachrichtigung vom Python-Server an Android
Herstellen einer Verbindung von Python zu MySQL unter CentOS 6.4
Portieren und Ändern des Doublet-Solvers von Python2 auf Python3.
Zugriff auf RDS von Lambda (Python)
[Python] Fluidsimulation: Implementieren Sie die Diffusionsgleichung
Python> Ausgaben von 1 bis 100, 501 bis 600> Für CSV
Einführung in die diskrete Ereignissimulation mit Python # 2
In Python von Markdown in HTML konvertieren
[Amazon Linux] Wechsel von der Python 2-Serie zur Python 3-Serie
API-Erklärung zum Berühren von Mastodon aus Python
Einführung in Vektoren: Lineare Algebra in Python <1>
Stellen Sie von Python aus eine Verbindung zur Websocket-API von coincheck her
Senden Sie eine Nachricht von Slack an einen Python-Server
Bearbeiten Sie Excel in Python, um eine Pivot-Tabelle zu erstellen
Auf Python 2.7.9 aktualisiert
So öffnen Sie einen Webbrowser über Python
Einführung in OPTIMIZER ~ Von der linearen Regression über Adam bis Eva
Studie aus Python Hour7: Verwendung von Klassen
[Python] Konvertieren von DICOM in PNG oder CSV
Excel-Datei aus Python importieren (in DB registriert)
Ich möchte mit Python eine E-Mail von Google Mail senden.
[Python] Ich möchte 7DaysToDie von Discord aus verwalten! 1/3
Von der Datei zur Diagrammzeichnung in Python. Grundstufe Grundstufe
[Python] Lesen von Daten aus CIFAR-10 und CIFAR-100
[Python] Erstellen Sie eine Tabelle von Pandas DataFrame zu Postgres
[Bash] Holen Sie sich die Kraft von Python aus Bash mithilfe der folgenden Dokumentation
So generieren Sie ein Python-Objekt aus JSON
Summe von 1 bis 10
SQL zu SQL
Wie man gut mit Linux-Befehlen aus Python umgeht
Ambisonics Simulation Python
[Python] Fluss vom Web-Scraping zur Datenanalyse
MeCab von Python