[Python] Simulation de fluide: de linéaire à non linéaire

introduction

Je voudrais résumer les connaissances nécessaires à la construction d'un code d'analyse numérique des fluides pour les liquides, tout en étudiant également la dynamique numérique des fluides (CFD). Je pense que la dynamique des fluides numérique est une étude assez difficile, alors j'aimerais l'écrire d'une manière facile à comprendre pour les débutants. Il semble qu'il y ait beaucoup d'erreurs, alors veuillez nous contacter si vous le trouvez. Je vous serais également reconnaissant de bien vouloir commenter ce qui est difficile à comprendre. Nous le mettrons à jour de temps en temps.

Public cible

séries

Contenu approximatif de cet article

Comme étape préliminaire pour construire l'équation de base du fluide nécessaire à la simulation de l'eau, nous expliquerons comment gérer les équations non linéaires. Enfin, nous effectuerons une simulation numérique de l'équation de Bergers, qui est une équation non linéaire.

Bibliothèque utilisée: Numpy, Scipy, Matplotlib

Je vais faire cette version bidimensionnelle avec un tel gars.

burgers.gif

table des matières

chapitre Titre Remarques
1. Linéaire et non linéaire
1-2. Caractéristiques des équations linéaires et non linéaires
1-3. Comment faire la distinction entre linéaire et non linéaire Bean bien informé
2. Implémentation de l'équation de diffusion par transfert linéaire
2-1. Équation de diffusion par transfert linéaire
2-2. la mise en oeuvre
3. Implémentation de l'équation de Bergers
3-1. Quelle est l'équation de Burgers?? Correspondance avec l'équation de diffusion de transfert non linéaire
3-2. la mise en oeuvre
3-3. Complément à l'expression discrète
4. Extension à 2 dimensions
4-1. Equation de transfert bidimensionnelle
4-2. Équation de diffusion bidimensionnelle
4-3. Equation bidimensionnelle de Burgers

1. Linéaire et non linéaire

1-1. Caractéristiques des équations linéaires et non linéaires

Ce qui suit est un bref résumé des équations linéaires et non linéaires.

équation Fonctionnalité
linéaire 1.Une relation proportionnelle tient
2. Peut trouver une solution
non linéaire 1.La relation proportionnelle ne tient pas
2. En gros, je ne trouve pas de solution(Seules quelques équations telles que l'équation kdV peuvent être résolues)

Les équations fluides à traiter à l'avenir sont classées comme des équations non linéaires, et fondamentalement aucune solution ne peut être obtenue. Par conséquent, lors de la recherche de la solution d'une équation non linéaire, la politique de base est d'effectuer un calcul numérique en utilisant les équations discrètes introduites dans les articles précédents et de trouver la solution approximativement. ..

1-2. Comment faire la distinction entre linéaire et non linéaire

En dehors du sujet, je voudrais expliquer comment distinguer les équations linéaires et non linéaires. Pour faire la distinction entre les équations linéaires et non linéaires, déterminez si la superposition de solution tient.

Je voudrais expliquer en utilisant l'équation de translocation présentée dans Article précédent comme exemple.

L'équation Advance est $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $ Il était représenté par. Cependant, le taux de transfert c est une constante (ce qui signifie que la substance se déplace à une vitesse constante de c).

Lorsqu'il y a deux solutions $ u_1 $ et $ u_2 $ dans cette équation, respectivement

\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

Est établi. Par conséquent, pour $ u_1 + u_2 , qui est la somme des deux solutions, $ \frac{\partial (u_1+u_2)}{\partial t} + c \frac{\partial (u_1+u_2)}{\partial x} = 0 $$ Est établie, et nous pouvons voir que cette équation est linéaire.

Maintenant, changeons la vitesse de transfert constante c en la 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

Alors, dans cette équation, $ u_1 + u_2 $, qui est la somme des deux solutions, n'est pas une solution, nous savons donc que c'est une équation non linéaire.

Dans cet article, je voudrais traiter de telles équations non linéaires.

2. Mise en œuvre de l'équation de diffusion par transfert linéaire

2-1. Équation de diffusion par transfert linéaire

Tout d'abord, je voudrais traiter de l'équation de diffusion par transfert linéaire, qui sert également de signification développementale à la série jusqu'à présent. Cette équation est comme une combinaison d'une équation de translocation (une équation qui signifie le mouvement d'une substance) et une équation de diffusion (une équation qui signifie l'uniformité des quantités physiques) comme le montre l'équation ci-dessous (c et). $ \ Nu $ est une constante). À propos, la version non linéaire de ceci est l'équation de Burgers qui apparaîtra dans le prochain chapitre.

2-2. Mise en œuvre

Pour le moment, nous allons l'implémenter en utilisant la méthode CIP et la méthode itérative non stationnaire. En termes simples, la méthode CIP est une méthode de prédiction des valeurs futures à partir de la valeur actuelle et de sa valeur différentielle. Pour plus de détails, voir l 'article et ce pdf Veuillez consulter jp / netlab / mhd_summer2001 / cip-intro.pdf). La méthode itérative non stationnaire est une méthode pour trouver la solution d'équations simultanées unidimensionnelles en changeant les coefficients. Si vous voulez en savoir plus, j'ai écrit Article sur les équations de diffusion. Voir / items / 97daa509991bb9777a4a) et Article sur la bibliothèque d'itération non stationnaire de Python.

Ce n'est pas si difficile à faire. Dans la méthode CIP, la phase de translocation et la phase de non-translocation sont considérées séparément. En bref, il vous suffit de calculer u à partir de l'équation de translocation et d'utiliser le résultat du calcul u pour résoudre l'équation de diffusion. En termes d'image, la vitesse de déplacement par translocation et la vitesse d'uniformisation par diffusion? (Fréquence?) Sont différentes, donc je pense que c'est comme penser séparément. L'équation de diffusion de translocation peut être exprimée comme suit en la divisant en une phase de translocation et une phase de non-translocation. Lors de l'utilisation de la méthode CIP, des informations sur la valeur différentielle $ \ frac {\ partial u} {\ partial x} $ sont également requises, elles sont donc effectuées en même temps que le calcul de u.

\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  #La valeur limite a été définie sur 0
    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
#Développer le temps
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)  #La valeur limite a été définie sur 0
    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. Mise en œuvre de l'équation de Bergers

3-1 Qu'est-ce que l'équation de Burgers?

C'est l'une des équations aux dérivées partielles non linéaires, et c'est une version non linéaire de l'équation de diffusion de translocation expliquée dans le chapitre précédent. Si vous faites des hypothèses sur l'équation fluide et que vous la simplifiez considérablement, cela devient l'équation de Burgers. Cette équation est une équation non linéaire, mais il est connu qu'une solution exacte peut être obtenue en utilisant une transformation de saut d'appel. Par conséquent, il est souvent utilisé pour examiner la validité des méthodes de calcul numérique.

3-2. Mise en œuvre

Implémente l'équation Burgers. Ceci est également implémenté par la méthode CIP. Elle est la suivante lorsqu'elle est divisée en une phase de translocation et une phase de non-translocation.

\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.

Il n'y a pas beaucoup de changements par rapport à l'équation de diffusion par transfert linéaire. Où changer

J'y pense. Nous avons apporté quelques modifications à la dispersion, voir la section suivante pour plus de détails.

Si vous considérez brièvement les résultats du calcul, vous pouvez voir que plus la hauteur est élevée, plus la vitesse est rapide, de sorte que le côté supérieur de la vague rattrape le côté inférieur sur le chemin, formant une surface de discontinuité abrupte semblable à une onde de choc. Après cela, vous pouvez voir que la solution devient progressivement terne en raison de l'effet de diffusion.

burgers.gif

def solve_advection_cip(u_cip, partial_u_cip, lower_boundary, dt=Delta_t, velocity=C):
    """
Résolution de l'équation de transfert par la méthode CIP
    """
    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):
    """
Résolvez l'équation de diffusion.
    
    Parameters
    ----------
    u_array : array-like
vecteur n
    upper_boundary : numpy.float64
        u_array[n]À côté de
    lower_boundary : numpy.float64
        u_array[0]À côté de
    diffusion_number : numpy.float64
Numéro de propagation. Nombre positif.
    update_partial : Bool, default False
Défini sur True lors de la mise à jour de la formule différentielle.
    u_array_integral : array-like, default None
Utilisé lors de la mise à jour de la formule différentielle.[lower_boundary, u_array, upper_boudary]N sous la forme de+Vecteur de 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Méthode semi-implicite. Manivelle-Nicolson scheme
        1: Fully implicit method.Différence de temps de retraite.

    Returns
    -------
    u_array : array-like
matrice à n lignes
    """
    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:
        #Faites cela lors de la mise à jour du différentiel
        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
#Définition des conditions de stabilité
cfl_condition = 1
diffusion_number_restriction = 1/2
#Développer le temps
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:
        #Jugement des conditions de la LCF
        # cfl_S'il est plus grand que la condition, réduisez la largeur du pas de temps dt afin que la condition CFL soit satisfaite.
        delta_t = cfl_condition * Delta_x / np.max(np.abs(u_cip))
    if diffusion_number > diffusion_number_restriction:
        #Jugement du nombre de diffusion
        # diffusion_number_S'il est plus grand que la restriction, réduisez la largeur du pas de temps dt afin que la condition du nombre de diffusion soit satisfaite.
        delta_t = diffusion_max * Delta_x ** 2 / Nu
    #Résoudre l'équation de translocation
    u_cip_star, partial_u_cip_star = solve_advection_cip(u_cip, partial_u_cip, u_lower_boundary, delta_t, u_cip[1:])
    #Résoudre l'équation de diffusion
    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)  #La valeur limite a été définie sur 0
    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. Formule de discrétion de la phase de diffusion

Équation de diffusion

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

Est exprimé par une expression discrète

\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)

Ce sera. Considérant l'avenir, contrairement à Article sur l'équation de diffusion, l'expression discrète est exprimée en utilisant $ \ lambda $. Pour résumer la valeur de ce $ \ lambda $

\lambda Méthode de discrimination Fonctionnalité
0 Différence de centre Méthode explicite. Calculez le décalage horaire en utilisant uniquement la valeur de l'heure actuelle.
1/2 Méthode implicite de Crank Nicholson Méthode semi-implicite. Calculez la différence de temps en utilisant la moitié de la valeur actuelle et la moitié de la valeur la fois suivante.
1 Différence de temps de retraite Méthode complètement implicite. Calculez le décalage horaire en utilisant uniquement la valeur de temps suivante.

Ce sera. L'essentiel est que vous pouvez utiliser $ \ lambda $ pour implémenter trois techniques de dispersion différentes.

Si vous transformez cela et amenez la valeur du temps n + 1 sur le côté gauche $ -\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} $

En supposant que les points de grille existent dans la plage de 1 à M, la valeur de la limite est représentée par $ u_0, u_ {M + 1} $, et elle est convertie en un affichage matriciel.

  \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)

Si vous voulez mettre à jour $ \ frac {\ partial u} {\ partial x} $, changez u dans cette matrice en $ \ frac {\ partial u} {\ partial x} $ et déplacez-le vers la droite.

\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)

Ajouter.

Ceci est implémenté dans la section précédente. La méthode BiCG STAB a été utilisée comme méthode de solution itérative. La raison en est que le nom est cool.

4. Extension à 2 dimensions

Je veux dessiner une image sympa, donc je la ferai en deux dimensions dans ce chapitre. C'est juste une expansion de formule.

4-1. Équation de transfert bidimensionnelle

Comme dans l'exemple, nous allons l'implémenter en utilisant la méthode CIP. Pour une extension détaillée des formules, reportez-vous à la Méthode CIP 2D dans la référence. Puisqu'il y a 10 inconnues, il suffit de créer 10 équations. Vous n'avez pas à penser aux coordonnées (i + 1, j + 1) au moment de la différenciation. Les références font un bon travail pour réduire la quantité de calcul, mais c'est déroutant, donc je vais le faire honnêtement ici.

Tout comme dans une 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

Définissez $ F (X, Y) $ avec complétion cubique comme dans. À partir de la condition que ** la valeur de la fonction et la valeur différentielle sont continues sur les points de grille **

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}

Les 10 formules de Cependant, la valeur au point de grille (i, j) de la fonction f et la valeur différentielle sont respectivement exprimées en $ f_ {i, j}, \ frac {\ partial f_ {i, j}} {\ partial x} $. Je vais. En remplaçant ceci dans la formule de $ F (X, Y) $ pour trouver le coefficient, $ 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} $

De ce qui précède, lorsque la vitesse est constante, $ \ frac {\ partial u} {\ partial t} + c_1 \ frac {\ partial u} {\ partial x} + c_2 \ frac {\ partial u} {\ partial y} = 0 $ satisfait également la valeur différentielle, donc la valeur et la valeur différentielle après la largeur du pas de temps $ \ Delta t (n + 1 pas) $ secondes sont

    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)

Si la vitesse n'est pas constante $ \ frac {\ partial u} {\ partial t} + u \ frac {\ partial u} {\ partial x} + v \ frac {\ partial u} {\ partial y} = 0 $ Le terme supplémentaire qui apparaît après la différenciation spatiale sera calculé comme un terme de non-translocation.

2d-advection.gif

4-2. Équation de diffusion bidimensionnelle

Équation de diffusion bidimensionnelle

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

Est exprimé par une expression discrète

  \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\}

Triez les grilles dans une ligne pour faciliter cette implémentation dans votre code. Définissez le point de grille de l'ensemble du système sur $ M = NX \ times NY $, et définissez la position du point de grille (i, j) sur le mth (début 0 selon la notation Python) à partir du point de grille (0,0). Je vais. Ensuite, le point de grille (i + 1, j) est le m + 1, et le point de grille (i, j + 1) est le m + NXth. En appliquant cela à la formule ci-dessus, si vous amenez la valeur du temps n + 1 sur le côté gauche

- 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}

Donc, si vous l'implémentez en incluant la valeur limite, vous avez terminé. Comme l'affichage matriciel est long, je vais l'ignorer ici.

2d-diffusion.gif

4-3. Équation bidimensionnelle de Burgers

Rendons-le bidimensionnel pour le moment. Changez la notation et changez u en f pour faire une distinction avec la vitesse.

\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.

Mettez simplement cela en œuvre. Comme indiqué ci-dessous, je ne sais pas si cela correspond, mais vous pouvez obtenir un résultat de calcul comme celui-ci.

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):
    """
Résolvez l'équation de diffusion. Il comprend également des éléments de termes de non-translocation.
    
    Parameters
    ----------
    u_2d : array-like
matrice de lignes nx × ny
    upper_boundary : numpy.float64
        u_array[n]À côté de
    lower_boundary : numpy.float64
        u_array[0]À côté de
    d_number_* : numpy.float64
Numéro de propagation. Nombre positif.
    update_partial : Bool, default False
Défini sur True lors de la mise à jour de la formule différentielle.
    u_array_integral : array-like, default None
Utilisé lors de la mise à jour de la formule différentielle.[lower_boundary, u_array, upper_boudary]N sous la forme de+Vecteur de 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Méthode semi-implicite. Manivelle-Nicolson scheme
        1: Fully implicit method.Différence de temps de retraite.

    Returns
    -------
    u_2d : array-like
matrice de lignes nx × ny
    """
    #Obtenir la taille de la matrice
    nx, ny = u_2d.shape
    matrix_size = nx * ny
    # Ax=Créer A et b pour 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)
    #Définir les conditions aux limites/ 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:
        #Exécuté lors de la mise à jour de l'expression différentielle dans la direction x. Non utilisé lors de la résolution d'équations de diffusion ordinaires.
        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:
        #Exécuté lors de la mise à jour de l'expression différentielle dans la direction y. Non utilisé lors de la résolution d'équations de diffusion ordinaires.
        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)

#Équation de Bugers 2D
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)

#Définition des conditions de stabilité
cfl_condition = 1
diffusion_number_restriction = 1/2
#Développer le temps
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:
        #Jugement des conditions de la LCF
        # cfl_S'il est plus grand que la condition, réduisez la largeur du pas de temps dt afin que la condition CFL soit satisfaite.
        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:
        #Jugement du nombre de diffusion
        # diffusion_number_S'il est plus grand que la restriction, réduisez la largeur du pas de temps dt afin que la condition du nombre de diffusion soit satisfaite.
        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) #Pour le moment
    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())

Résumé

La prochaine fois, je voudrais traiter des équations de base pour les fluides. De plus, le calcul est devenu assez lourd, je voudrais donc l'accélérer si j'en ai envie.

Les références

Recommended Posts

[Python] Simulation de fluide: de linéaire à non linéaire
Changements de Python 3.0 à Python 3.5
Changements de Python 2 à Python 3.0
Publier de Python vers Slack
Flirter de PHP à Python
Anaconda mis à jour de 4.2.0 à 4.3.0 (python3.5 mis à jour vers python3.6)
Passer de python2.7 à python3.6 (centos7)
Connectez-vous à sqlite depuis python
Appelez Matlab depuis Python pour optimiser
Publication de Python sur la chronologie Facebook
[Lambda] [Python] Publier sur Twitter depuis Lambda!
Connectez-vous à la base de données utf8mb4 à partir de python
Python (de la première fois à l'exécution)
Publier une image de Python sur Tumblr
[Python] Simulation de fluide: équation de Navier-Stokes incompressible
Comment accéder à wikipedia depuis python
Python pour passer d'une autre langue
N'a pas changé de Python 2 à 3
Mettre à jour Mac Python de 2 à 3
De Python à l'utilisation de MeCab (et CaboCha)
Introduction à la simulation d'événements discrets à l'aide de Python # 1
[Python] Simulation de fluides: implémenter l'équation de transfert
Comment mettre à jour Google Sheets à partir de Python
Manuel Python privé (mis à jour de temps en temps)
Je veux utiliser jar de python
Conversion de katakana en voyelle kana [python]
Notification push du serveur Python vers Android
Connexion de python à MySQL sur CentOS 6.4
Portage et modification du solveur de doublets de python2 vers python3.
Comment accéder à RDS depuis Lambda (python)
[Python] Simulation de fluides: implémenter l'équation de diffusion
Python> Numéros de sortie de 1 à 100, 501 à 600> Pour csv
Introduction à la simulation d'événements discrets à l'aide de Python # 2
Convertir de Markdown en HTML en Python
[Amazon Linux] Passage de la série Python 2 à la série Python 3
Explication API pour toucher mastodonte de python
Introduction aux vecteurs: Algèbre linéaire en Python <1>
Connectez-vous à l'API Websocket de Coincheck depuis Python
Envoyer un message de Slack à un serveur Python
Modifier Excel à partir de Python pour créer un tableau croisé dynamique
Mis à jour vers Python 2.7.9
Comment ouvrir un navigateur Web à partir de python
Introduction à OPTIMIZER ~ De la régression linéaire à Adam à Eve
Étude de Python Hour7: Comment utiliser les classes
[Python] Conversion de DICOM en PNG ou CSV
Importer un fichier Excel depuis Python (enregistré dans DB)
Je souhaite envoyer un e-mail depuis Gmail en utilisant Python.
[Python] Je veux gérer 7DaysToDie depuis Discord! 1/3
Du dessin de fichier au graphique en Python. Élémentaire élémentaire
[Python] Comment lire les données de CIFAR-10 et CIFAR-100
[python] Créer une table de pandas DataFrame vers postgres
[Bash] Obtenez la puissance de python de bash en utilisant la documentation ici
Comment générer un objet Python à partir de JSON
Somme de 1 à 10
sql à sql
Comment bien gérer les commandes Linux à partir de Python
simulation ambisonique Python
[Python] Flux du scraping Web à l'analyse des données
MeCab de Python