[Python] Simulation de fluides: implémenter l'équation de transfert

introduction

Je voudrais résumer (en plusieurs articles) les connaissances nécessaires pour construire un code d'analyse numérique des fluides pour l'eau, tout en étudiant également la dynamique numérique des fluides (CFD), qui est une étude liée à la simulation de fluides tels que l'air et l'eau.

Je voudrais l'écrire pour que même les débutants puissent facilement le comprendre. 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 traiter les équations de base du fluide requises pour la simulation de l'eau, nous résumerons brièvement et implémenterons les équations de transfert.

Journal des modifications

Equation de translocation (formule représentant le mouvement des substances)

Nous résoudrons l '[équation de transfert] unidimensionnelle (http://zellij.hatenablog.com/entry/20140805/p1) par la méthode des différences finies.

Équation avancée $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $

Ce que signifie cette équation, c'est qu'une fonction u se déplace à la vitesse c. En bref, c'est une équation qui exprime que la quantité physique se déplace le long du flux de vitesse c.

advection_example.png

Dispersez ceci et résolvez l'équation. Cette fois, nous allons le résoudre par une méthode appelée Positive Solution. La méthode explicite est simplement une méthode de recherche de valeurs futures en utilisant uniquement les valeurs de temps actuelles. En plus de la méthode explicite, il existe également une méthode appelée méthode implicite, qui est une méthode de calcul basée sur la valeur actuelle et la valeur future.

La méthode explicite à mettre en œuvre cette fois est la suivante. J'ai choisi les quatre suivants à volonté.

Pour chacun, j'écrirai les expressions discrètes. Dans la formule, l'indice j signifie le lieu et n le temps. Par conséquent, $ f_j ^ n $ signifie la valeur de la fonction f située à l'emplacement $ x_j $ au moment $ n $.

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

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

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

  4. CIP(Constrained Interpolation Profile scheme) method

Les détails seront décrits plus tard. Vous pouvez également vous reporter à Introduction à la loi CIP.

Conditions stables de la méthode explicite

Je voudrais mentionner brièvement les conditions de stabilité importantes lors de l'utilisation de la méthode explicite. Numériquement, la vitesse de propagation de l'information $ \ frac {\ Delta t} {\ Delta x} $ (la vitesse obtenue en divisant la largeur de la grille par le temps) doit être supérieure ou égale à la vitesse physique du signal $ c . Sinon, le phénomène physique se déplacera plus rapidement que le calcul numérique et le phénomène physique ne pourra pas être exprimé numériquement. Basé sur cette idée $ c \leq \frac{\Delta x}{\Delta t} $ La condition est dérivée. Transformé cela $ \nu = c \frac{\Delta t}{\Delta x} \leq 1 $$ Est appelée la condition CFL et est utilisée comme condition numériquement stable. Le côté gauche est également appelé le numéro Coolan. En bref, la largeur du pas de temps $ \ Delta t $ doit être inférieure ou égale à $ \ frac {\ Delta x} {c} $. Notez que ces conditions CFL ne sont pas pertinentes lors de l'utilisation de la méthode implicite, donc la méthode implicite a l'avantage de pouvoir augmenter les incréments de temps par rapport à la méthode explicite. Cependant, le calcul se complique en conséquence.

la mise en oeuvre

Le langage utilisé est Python. Je vais l'écrire sans utiliser de nombreuses fonctions Numpy pour que la formule de calcul soit facile à comprendre.

Objectif de calcul

Calculez le transfert d'une onde rectangulaire. La largeur de la grille $ \ Delta x $, la vitesse physique du signal $ c $ et la largeur du pas de temps $ \ Delta t $ sont toutes calculées comme des constantes. Cette fois, il est calculé comme $ \ Delta x = 1 $, $ c = 1 $, $ \ Delta t = 0,2 $ selon la référence. Par conséquent, il est fixé à $ CFL = 0,2 $ et la condition CFL est satisfaite.

advection_answer.png 
# Creat square wave
Num_stencil_x = 101
x_array = np.arange(Num_stencil_x)
u_array = np.where(((x_array >= 30) |(x_array < 10)), 0.0, 1.0)
u_lower_boundary = 0.0
u_upper_boundary = 0.0
Time_step = 200
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
total_movement = C * Delta_t * (Time_step+1)
exact_u_array = np.where(((x_array >= 30 + total_movement) |(x_array < 10 + total_movement)), 0.0, 1.0)
plt.plot(x_array, u_array, label="Initial condition")
plt.plot(x_array, exact_u_array, label="Answer")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)
  1. FTCS(Forward in Time and Central Difference in Space) $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) $

answer_ftcs.png

La solution vibre. Je n'ai pas du tout de réponse

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

2. Différenciation au vent

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

answer_upwing.png

La diffusion numérique est grande et la solution devient lisse.

u_upwind = u_array.copy()
#Définir le pas de temps
for n in range(Time_step):
    u_old = u_upwind.copy()
    u_upwind[0:1] = u_old[0] - CFL * (u_old[1] - u_lower_boundary)
    for j in range(1,Num_stencil_x):
        u_upwind[j:j+1] = u_old[j] - CFL * (u_old[j] - u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_upwind, label="Upwind differencing")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Upwind differencing)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)
  1. Lax-Wendroff scheme $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) + \frac{c^2}{2} \left( \frac{\Delta t}{\Delta x} \right)^2 \left(u_{j+1}^n - 2 u_j^n+ u_{j-1}^n \right) $

answer_laxw.png

Il y a une vibration numérique et une solution terne.

u_lw= u_array.copy()
#Définir le pas de temps
for n in range(Time_step):
    u_old = u_lw.copy()
    u_lw[0:1] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary) \
                    + CFL**2 / 2 * (u_old[1] - 2 * u_old[0] + u_lower_boundary)
    u_lw[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1]) \
                    + CFL**2 / 2 * (u_upper_boundary - 2 * u_old[-1] + u_old[-2])
    for j in range(1,Num_stencil_x-1, 1):
        u_lw[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1]) \
                    + CFL**2 / 2 * (u_old[j+1] - 2 * u_old[j] + u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_lw, label="Lax-Wendroff scheme")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Lax-Wendroff scheme)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)
  1. CIP(Constrained Interpolation Profile scheme) method

L'idée de base de la méthode CIP est que lorsqu'une fonction f est d'abord discrétisée, les quantités physiques d'un point $ x_j $ et d'un point $ x_ {j-1} $ à côté d'elle sont reliées en douceur. En faisant cela, les valeurs différentielles des grandeurs physiques respectives $ f_j $ et $ f_ {j-1} $ peuvent être utilisées pour l'évolution temporelle sans diverger l'inclinaison entre les réseaux.

En d'autres termes, à condition que la ** valeur de la fonction et la valeur différentielle soient continues sur les points de la grille **, la relation entre les deux points de la grille peut être exprimée par la fonction complémentaire cubique $ F (x) $.

Dans la zone de $ x_ {j-1} <= x <x_j $ $ F(x) = a X^3 + b X^2 +c X + d \quad ,where \quad X = x - x_j $ Si $ F (x) $ est défini par complétion cubique comme dans, ** à partir de la condition que la valeur de la fonction et la valeur différentielle soient continues sur le point de grille **,

F(0) = f_j^n, \quad F(\Delta x) = f_{j+1}^n, \quad \frac{\partial F(0)}{\partial x} = \frac{\partial f_j^n}{\partial x}, \quad \frac{\partial F(\Delta x)}{\partial x} = \frac{\partial f_{j+1}^n}{\partial x} \\ , where \quad \Delta x = x_{j+1} - x_j

Est satisfait, donc si vous remplacez ceci dans la formule de $ F (x) $ et trouvez $ a, b, c, d $

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

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

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

answer_cip.png

La diffusion numérique et les vibrations sont extrêmement plus petites que les trois autres méthodes. ~~ Cependant, le transfert n'était pas aussi bon que Références. ~~ (2020/02/20 Après avoir apporté les corrections que j'ai reçues dans les commentaires, cela s'est considérablement amélioré.)

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_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = 0
    partial_u_cip[0] = 0
    for j in range(1, Num_stencil_x):
        a = (partial_u_old[j] + partial_u_old[j-1]) / Delta_x**2 - 2.0 * (u_old[j] - u_old[j-1]) / Delta_x**3
        b = 3 * (u_old[j-1] - u_cip[j]) / Delta_x**2 + (2.0*partial_u_old[j] + partial_u_old[j-1]) / Delta_x
        c = partial_u_old[j]
        d = u_old[j]
        xi = - C * Delta_t  # C > 0
        u_cip[j:j+1] = a * xi**3 + b * xi**2 + c * xi + d
        partial_u_cip[j:j+1] = 3 * a * xi**2 + 2 * b * xi + c
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_cip, label="CIP method")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(CIP)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

Résumé

Les bases des équations de translocation qui composent l'équation de fluide sont résumées avec une brève mise en œuvre. Les quatre suivants ont été mis en œuvre.

La prochaine fois, je parlerai de Équation de diffusion.

Les références

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

Recommended Posts

[Python] Simulation de fluides: implémenter l'équation de transfert
[Python] Simulation de fluides: implémenter l'équation de diffusion
[Python] Simulation de fluide: équation de Navier-Stokes incompressible
Implémenter le modèle Singleton en Python
[Python] Observons le vortex de Kalman par la méthode de Boltzmann en treillis. [Simulation de fluide]
simulation ambisonique Python
Trouvez la solution de l'équation d'ordre n avec python
Résolution d'équations de mouvement en Python (odeint)
Mettre en œuvre REPL
Implémenter la solution de l'algèbre de Riccati en Python
Mettre en œuvre des recommandations en Python
le zen de Python
Implémenter XENO avec python
Implémenter sum en Python
[Python] Fractionner la date
Implémenter Traceroute dans Python 3
[Python] Régression LASSO avec contrainte d'équation utilisant la méthode du multiplicateur
Avoir le graphique d'équation de la fonction linéaire dessiné en Python
J'ai essayé d'implémenter la fonction d'envoi de courrier en Python