[Python] Simulation de fluide: équation de Navier-Stokes incompressible

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). Il semble qu'il y ait beaucoup d'erreurs, alors veuillez nous contacter si vous le trouvez.

Public cible

séries

Contenu approximatif de cet article

Cet article traite des simulations numériques monophasées telles que gaz uniquement ou liquide uniquement. Après avoir brièvement expliqué l'équation gouvernante du fluide, il est configuré pour mettre en œuvre l'algorithme de calcul numérique du fluide incompressible nécessaire à la simulation du liquide. Enfin, nous simulerons l'écoulement de la cavité (sensation de rotation du fluide dans la zone / figure ci-dessous). En ce qui concerne le flux de cavité, @ eigs Solution numérique de l'équation incompressible de Navier-Stokes 1: Introduction est facile à comprendre et le graphique est également beau.

2d-cavity.gif

table des matières

chapitre Titre Remarques
1. Solution fluide
1.1. Equation de régulation fluide
1.2. Compressible et incompressible
2. Algorithme de calcul de fluide incompressible
2.1. Méthode MAC Marker And Cell method
2.2. Méthode par étapes fractionnaires Aussi appelée méthode de l'étape partielle
2.3. Méthode SIMPLE Semi-Implicit Method for Pressure-Linked Equation method
2.4. Méthode CCUP CIP-Combined Unified Procedure method
2.5. Différence spatiale Concernant le placement variable
3. la mise en oeuvre Méthode CCUP

1. Solution fluide

1-1. Équation de contrôle du fluide

Un fluide se compose de trois équations: une équation continue, une équation Navigator Stokes et une équation d'énergie. Chaque équation est dérivée de la conservation de masse, de la conservation de l'élan et de la conservation de l'énergie. C'est peut-être plus déroutant, mais j'ajouterai également une image de chaque formule.

a. Expression continue

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

Il peut être dérivé du fait que le changement de masse dans un certain volume de test est égal au débit massique $ \ rho \ vec {u} $ qui entre et sort de la surface divisé par l'aire. Pour expliquer avec un exemple facile à comprendre, la formule ci-dessus montre que le montant de la monnaie dans le portefeuille est égal au montant obtenu en soustrayant le montant d'argent passé par la porte du portefeuille.

money.png

b. Équation de Navier Stokes (équation de mouvement fluide)

\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} \\
 \quad or \\
\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \nabla \cdot {\bf T}_{\nu} + \vec{F_B} \\
, where \quad {\bf T}_{\nu} = \lambda (\nabla \cdot \vec{u}) {\bf I} + 2 \mu {\bf D} \quad :Tenseur des contraintes visqueuses\\
{\bf D}  = \frac{1}{2} \left\{(\nabla \vec{u}) + (\nabla \vec{u})^T \right\} :Tenseur du taux de déformation\\
\quad \vec{F_B} \quad :Force de volume(Principalement la gravité) \\
\quad \mu :Coefficient visqueux

La quantité de mouvement du fluide est $ gradient de pression: \ nabla p $, $ diffusion due à la viscosité: \ nabla \ cdot {\ bf T} _ {\ nu} $, $ force volumique (gravité etc.): \ vec {F_B} $ Indique que vous serez affecté. Pour donner un exemple d'image, il peut être utile d'imaginer que lorsqu'un garçon frappe une balle, la vitesse de la balle change en raison des effets de la différence de pression, de la viscosité et de la gravité.

navier_stokes.png

c. Équation énergétique

 \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L} \\
, where \quad  \dot{L} :Quantité de chaleur générée par les réactions chimiques, etc.\\
\quad \kappa :Conductivité thermique

cependant,

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

La quantité de changement de l'énergie interne du fluide est le travail $ \ nabla \ cdot \ vec {u} $ dû à la pression, le bilan énergétique dû à la conduction thermique $ \ nabla \ cdot (\ kappa \ nabla T) $, et l'énergie dissipée due à la viscosité $ {\ bf Cela signifie qu'il est affecté par $ \ dot {L} $ tel que T} _ {\ nu} \ cdot \ nabla \ vec {u} $ et des réactions chimiques. Cette formule vous permet de calculer le changement de température du fluide. En bref, le changement de température du fluide est égal à la quantité de chaleur ajoutée.

energy_equation.png

Résumé

La figure ci-dessus peut être résumée comme suit.

flow_equations.png

L'équation gouvernante du fluide dans un système non conservateur est illustrée ci-dessus, et si le côté droit est défini sur 0, cela signifie l'équation de transfert de densité, de vitesse et d'énergie interne (formule qui signifie le mouvement des substances). Un système non conservateur est une forme dans laquelle des variables existent en dehors de la différenciation, et il n'y a aucune garantie que les résultats du calcul satisfont à la loi de conservation, il est donc nécessaire de confirmer le caractère prudent lors de l'exécution de calculs numériques. Le contraire du système non conservateur est appelé le système conservateur, qui est une expression dans laquelle toutes les variables sont incluses dans le différentiel.

1-2. Compressible et non compressible

L'équation présentée dans la section 1-1 est l'équation dominante pour les fluides compressibles. Un fluide compressible est un fluide dont la vitesse d'écoulement est d'environ 0,3 fois la vitesse du son (le nombre de Mach est de 0,3), et un fluide dont l'écoulement est plus lent est généralement reconnu comme un fluide incompressible. En bref, si le fluide se déplace trop vite, le fluide lui-même se dilate et se contracte (changements de densité), le rendant compressible, et s'il est lent ou dense, le fluide ne se dilate pas ou ne se contracte pas et est traité comme incompressible. Par conséquent, le liquide est essentiellement traité comme non comprimé.

flow.png

La formule présentée dans la section 1-1 est transformée en une formule non compressée. Tout ce que nous avons à faire est de traiter la densité comme une constante plutôt qu'une variable pour simplifier la formule. Les équations de Navier Stokes et d'énergie sont simplifiées en utilisant des équations continues. De plus, nous ignorerons ici les effets des réactions chimiques.

a. Expression continue

\nabla \cdot \vec{u} = 0

b. Équation de Navier Stokes

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

c. Équation énergétique

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

Les équations des fluides compressifs et non compressibles peuvent être résumées comme suit.

Équation de base Compressibilité 非Compressibilité
a.Formule continue \frac{\partial \rho}{\partial t} + (\vec{u} \cdot \nabla) \rho = - \rho \nabla \cdot \vec{u} $ \nabla \cdot \vec{u} = 0$
b.Équation de Navier Stokes $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} $ $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + \vec{F_B} $
c.Équation énergétique $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L}$ $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = - \nabla \cdot (\kappa \nabla T) + 2 \mu {\bf D}:{\bf D}$

2. Algorithme de calcul des fluides incompressibles

Puisque nous cherchons à simuler des liquides, ce chapitre décrira comment résoudre les fluides incompressibles. Fondamentalement, nous allons construire une équation discrète en utilisant une équation incompressible. Pour éviter les complications, nous ignorerons ici les équations d'énergie et entendons résoudre "explicitement" avec la méthode explicite et "résoudre implicitement" avec la méthode implicite. En gros, je me réfère à Compressible Fluid Dynamics.

Les algorithmes de calcul typiques sont résumés dans le tableau ci-dessous. La méthode MAC (Maker And Cell), la méthode Fractional Step et la méthode SIMPLE sont les algorithmes de calcul de base pour les fluides incompressibles. Puisqu'il est implémenté dans divers sites et livres, je n'expliquerai que le principe ici. À la fin de ce chapitre, nous mentionnons la méthode CCUP (CIP-Combined Unified Procedure) et l'implémentons dans le chapitre suivant. La méthode CCUP est une méthode de résolution de l'équation des fluides par la méthode CIP utilisée dans les articles précédents.

algorithme Terme de translocation Terme de diffusion Terme de pression
Système MAC Méthode explicite Méthode explicite Méthode implicite
Méthode par étapes fractionnaires Méthode explicite Méthode implicite Méthode implicite
---- Méthode implicite Méthode explicite Méthode implicite
Système de méthode SIMPLE Méthode implicite Méthode implicite Méthode implicite
Méthode CCUP Méthode explicite Méthode expliciteorMéthode implicite Méthode implicite

2-1. Méthode MAC

La méthode MAC (Maker And Cell) est une solution de base dans le calcul des fluides car elle permet de résoudre différents flux de manière stable. La procédure de calcul est la suivante.

  1. Résolvez l'équation de pression par la méthode implicite pour trouver la pression $ p ^ {n + 1} $.
  2. Avancez l'équation de vitesse au fil du temps pour trouver la vitesse $ \ vec {u ^ {n + 1}} $.

L'équation de pression et l'équation de vitesse sont calculées comme suit. Multipliez l'équation de Navier Stokes par $ \ nabla $ pour prendre la divergence et la disperser.

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

Il peut être obtenu. De plus, si le terme de translocation et le terme de viscosité de l'équation de Navier Stokes sont explicitement séparés dans le temps,

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

Il peut être obtenu. En résolvant délibérément $ \ nabla \ cdot \ vec {u ^ n} $ sans le mettre à 0, nous concevons des moyens de supprimer les erreurs lors de la dispersion.

Il existe la méthode SMAC (Simplified MAC), la méthode HSMAC (Highly Simplified MAC), etc. comme algorithme de calcul du système de méthode MAC. Ceux-ci ont été proposés pour résoudre la méthode MAC à grande vitesse.

2-2. Méthode par étapes fractionnaires

Dans la méthode MAC, le terme visqueux est résolu par la méthode explicite, et il est grandement affecté par la limitation de largeur de pas de temps (condition CFL). La méthode de l'étape fractionnaire décrite dans cette section est une méthode qui résout le problème des conditions CFL en évaluant implicitement le terme de viscosité. Aussi connue sous le nom de méthode par étapes partielles. En tant que caractéristique, dans la formule liée à la vitesse utilisée dans la méthode MAC, la vitesse et la pression ont été mélangées sur le côté droit, mais dans la méthode Fractional Step, la vitesse et la pression sur le côté droit sont complètement séparées comme suit.

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

Comme procédure de calcul

  1. Trouvez la vitesse provisoire $ \ bar {u} $ à partir de la formule supérieure de la formule de vitesse.
  2. Utilisez l'équation de pression pour trouver la pression $ p ^ {n + 1} $.
  3. Trouvez la vitesse $ \ vec {u ^ {n + 1}} $ en utilisant la deuxième formule à partir du haut de la formule pour la vitesse.

2-3. Méthode SIMPLE

Abréviation de la méthode semi-implicite pour la méthode d'équation liée à la pression, qui est formulée en fonction de la méthode implicite. vouloir en savoir davantage! Bases de l'analyse thermo-fluide 64 Chapitre 6 Méthode d'analyse thermo-fluide: 6.5.4 Méthode SIMPLE explique en détail.

  1. Trouvez implicitement la vitesse provisoire $ \ bar {u} $.
\bar{u} = \vec{u^n} - \Delta t \left\{(\bar{u} \cdot \nabla)\bar{u} + \nabla p - \nu \nabla^2 \bar{u}\right\}
  1. Mettez à jour la pression à l'aide de l'équation de pression.
\nabla^2 \delta p = \frac{\nabla \cdot \bar{u}}{\Delta t}\\ p^{n+1} = p + \delta p
  1. Mettez à jour la vitesse $ \ vec {u} $ avec la pression calculée.
\vec{u^{n+1}} - \bar{u} = - \Delta t \nabla(p^{n+1} - p)

2-4. Méthode CCUP

La méthode CCUP est une méthode construite par Yabe et al., Qui ont conçu la méthode CIP, comme lorsque l'équation de Burgers dans article précédent a été résolue par la méthode CIP. , L'équation gouvernante du fluide est calculée séparément pour la phase de translocation et la phase de non-translocation. Shimizu et al. et Himeno et al. Le papier de S est facile à comprendre. Puisque l'ordre de résolution des termes de translocation et de non-translocation diffère d'une personne à l'autre, nous les résoudrons dans l'ordre présenté dans le dernier article ci-dessous. En tant qu'algorithme

  1. Dans la phase de transfert, transférez $ densité: \ rho $, $ vitesse: \ vec {u} $, $ énergie interne: e $, $ pression: p $.
  2. Passez à la phase de non-transfert.
  3. Mettre à jour la vitesse, l'énergie interne et la pression à l'étape de diffusion. (Densité constante)
  4. Résolvez l'équation de pression implicitement et mettez à jour la pression.
  5. Obtenez la densité, la vitesse et l'énergie interne en utilisant une nouvelle pression.

Je vais résoudre dans l'ordre. Bref, après le transfert et la diffusion de la substance, je pense que la vitesse de déplacement est extrêmement rapide (la même vitesse que la vitesse du son) et le Tsuji est ajusté.

ccup_method.png

2-5. Différence spatiale

Avant d'entrer dans l'implémentation, nous pouvons l'implémenter sans connaître la méthode de la différence spatiale, mais je la mentionnerai brièvement. En effet, une vibration de pression appelée instabilité en damier peut se produire en fonction de l'endroit où la variable est placée dans l'espace.

Treillis régulier

regular_stencil.png

Treillis échelonné

staggered_grid.png

Réseau Colocate

collocated_grid_2.png

Dans la méthode ci-dessus, la politique de base est d'utiliser la grille échelonnée pour le calcul, mais dans la mise en œuvre de cet article, le calcul numérique est effectué en utilisant la grille régulière (~~ Je ne pouvais pas comprendre la motivation lors de l'utilisation de la grille échelonnée, etc. Alors ~~).

3. Mise en œuvre

Il est implémenté par la méthode CCUP. Le problème à résoudre est l'écoulement de la cavité, qui donne l'impression de donner une vitesse latérale à la paroi supérieure, provoquant un vent dans le sens des aiguilles d'une montre dans la zone. Le code est publié sur Github. Ignorez l'équation énergétique. Puisque le terme de transfert et le terme de diffusion ont été décrits dans l'article précédent, le nouveau contenu du code est la partie d'équation de pression, mais le calcul du terme de diffusion et La structure est presque la même. Dans le calcul lié à la pression, certains éléments tels que le terme de transfert sont omis. De plus, les conditions aux limites sont assez grossières

cavity_flow.png

Si vous exprimez la vitesse avec une ligne rouge et la pression avec un contour bleu, vous pouvez voir qu'un graphique comme celui-ci est généré.

2d-cavity.gif

scipy.sparse.dia_matrix

Nous utilisons scipy.sparse.dia_matrix pour créer la matrice d'équation simultanée unidimensionnelle $ {\ bf A} $ $ {\ bf A} \ vec {x} = \ vec {b} $. Cette matrice $ {\ bf A} $ est appelée matrice creuse, et comme la plupart des éléments de la matrice sont composés de 0, elle est spécifiée par scipy.sparse.dia_matrix Méthode de stockage à compression diagonale (DIA). est efficace en termes de mémoire et ainsi de suite. D'autres formats de stockage de matrice clairsemée sont décrits dans [Python] Articles qui permettent le calcul de matrices clairsemées à grande vitesse.

Un exemple est présenté ci-dessous.

#Préparez un ou plusieurs tableaux de même longueur. Notez que vous avez besoin de la même longueur.
data = np.array([np.arange(10), np.arange(20,30)], np.arange(30,40))  
#Déterminez le décalage de la matrice diagonale en fonction de l'indice des données. 0:Matrice diagonale,Nombre positif:Face supérieure,Nombre négatif:Inférieur
#Notez que si l'offset est positif, le tableau sera rogné depuis le début et si l'offset est négatif, le tableau sera rogné depuis l'arrière.
offsets = np.array([0, 2, -3])  
a_matrix = scipy.sparse.dia_matrix((data, offsets), shape=(len(data[0]), len(data[0])))  #Définissez la forme de la matrice avec la forme.

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

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

4. Résumé

Ensuite, je voudrais parler de la simulation biphasée gaz-liquide (gaz et liquide).

Les références

Recommended Posts

[Python] Simulation de fluide: équation de Navier-Stokes incompressible
[Python] Simulation de fluides: implémenter l'équation de transfert
[Python] Simulation de fluides: implémenter l'équation de diffusion
[Python] Simulation de fluide: de linéaire à non linéaire
simulation ambisonique Python
Première simulation de cellule nerveuse avec NEURON + Python
simulation ambisonique (problème externe) Python
Simulation de diffraction des rayons X avec Python