Essayez de dessiner un diagramme de contour en utilisant matplotlib avec OpenFOAM

introduction

Cet article concerne les astuces de OpenFOAM, qui est un logiciel open source pour l'analyse des fluides. Dans OpenFOAM, le résultat du calcul est dessiné en utilisant ParaView comme post-traitement (visualisation), mais fondamentalement l'image dessinée est au format raster et est bidimensionnelle. Le contour de la section transversale ne peut pas être édité au format vectoriel. De plus, lorsque la section transversale / variable à visualiser est décidée, le travail de démarrage de ParaView et de dessin devient difficile. Par conséquent, cet article résume comment dessiner de manière semi-automatique un contour de coupe bidimensionnelle à l'aide de matplotlib, qui est une bibliothèque de visualisation de python, et le générer au format vectoriel.

Environnement

Le logiciel requis cette fois-ci est OpenFOAM et Python. La version d'OpenFOAM a récemment été mise à jour majeure vers 4.0 ou v1606, mais cet article se concentre sur 2.4.x. Pour Python, 2.7 est la cible. Numpy, matplotlib, vtk (python-vtk), tkinter (python-tk) sont nécessaires en tant que packages Python, donc si ce n'est pas le cas, veuillez les installer en plus.

Après avoir construit l'environnement, vérifiez le fonctionnement avec le terminal.

OpenFOAM

$ foamInstallationTest
$ mkdir -p $FOAM_RUN
$ run
$ cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity ./
$ cd cavity
$ blockMesh
$ icoFoam

foamInstallationTest vérifie les paramètres d'environnement de base de l'OpenFOAM installé. Dans les commandes suivantes, le flux de cavité, qui est le didacticiel de base d'OpenFOAM, est exécuté, et s'il n'y a pas d'erreurs jusqu'à présent, c'est OK.

Python

$ python
>>> import numpy
>>> import matplotlib
>>> import vtk
>>> import Tkinter
>>> quit()

Si chaque package (module) est installé, rien ne sera affiché et l'invite sera renvoyée normalement. ʻImportError: No module Si vous obtenez une erreur nommée xxxx`, réinsérez le package approprié.

Acquisition de données de section transversale bidimensionnelle (VTK)

Tout d'abord, à partir du résultat de l'analyse d'OpenFOAM, acquérez les données de coupe bidimensionnelle que vous souhaitez dessiner au format vtk. Vous pouvez lire les données OpenFOAM avec ParaView et les données de section en sortie avec Slice, mais le but de cette fois est de dessiner automatiquement des diagrammes de contour avec des commandes ou des scripts autant que possible, donc ParaView qui nécessite une opération GUI est Non utilisé. Ici, nous utiliserons l'utilitaire OpenFOAM sample (qui est intégré dans postProcess dans OpenFOAM 4.0 et versions ultérieures).

Créez des données de coupe à l'aide des résultats du didacticiel sur la cavité effectué lors de la vérification des paramètres d'environnement. Tout d'abord, pour utiliser «sample», copiez et éditez le modèle de «sampleDict» qui est le fichier de paramétrage de cette commande.

$ run
$ cd cavity
$ cp $FOAM_UTILITIES/postProcessing/sampling/sample/sampleDict ./system
$ 

Le fichier de modèle original décrit diverses méthodes d'échantillonnage, y compris les commentaires d'utilisation, mais cette fois nous n'échantillonnerons que le plan $ x $ - $ y $ de la cavité bidimensionnelle, donc simple sampleDict comme suit. À

setFormat raw;
surfaceFormat vtk;
formatOptions
{
}
interpolationScheme cellPoint;
fields
(
    p
    U
);
sets
(
);
surfaces
(
    interpolatedPlane
    {
        type            plane;
        coordinateSystem
        {
            origin      (0.05 0.05 0.005);
            coordinateRotation
            {
                type    axesRotation;
                e1      (1 0 0);
                e2      (0 1 0);
            }
        }
        basePoint       (0 0 0);
        normalVector    (0 0 1);
        interpolate     true;
    }
);

Puisque «sets» est utilisé lors de l'acquisition de données unidimensionnelles telles que des lignes droites, laissez ce champ vide cette fois. Il existe plusieurs manières de spécifier des «surfaces», mais cette fois nous utiliserons celle obtenue par interpolation.

Après avoir édité sampleDict, exécutez la commande sample pour générer des données de section bidimensionnelles. Seules les données de la dernière fois (0.5) sont échantillonnées avec l'option -latestTime.

$ sample -latestTime

Lorsque la commande est exécutée, un répertoire appelé postProcessing / surfaces sera créé dans le répertoire courant, et les données vtk à l'heure spécifiée seront sorties.

Préparation du script Python

Le format du fichier vtk créé par l'utilitaire sample ci-dessus est ce que l'on appelle le format hérité version 2.0. En Python, si vous pouvez lire des données au format vtk et transmettre les données à matplotlib de la bibliothèque de visualisation, vous pouvez dessiner librement. Vous pouvez écrire le script vous-même, mais Alexey Matveichev a écrit [Tracer des fichiers VTK avec matplotlib et python-vtk](http://matveichev.blogspot.jp/2014/03/plotting-vtk-files-with-matplotlib] -et.html) est écrit, utilisez donc le script publié. Le script publié dans l'article a une version vtk de python-vtk de la série 5.x, donc si vous créez un environnement avec la nouvelle série 6.x, il est lié dans l'article [script publié sur github](https: Veuillez utiliser //gist.github.com/mrklein/44a392a01fa3e0181972). Ici, nous allons procéder en supposant que le système 5.x est utilisé.

Copiez le script d'origine et placez-le dans le répertoire d'analyse OpenFOAM (nommons-le plotUx.py) et modifiez-le un peu.

plotUx.py


#!/usr/bin/env python

import sys  # new
import os   # new
import numpy as N
import vtk
import matplotlib.pyplot as P

def load_velocity(filename):
    if not os.path.exists(filename): return None
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.ReadAllVectorsOn()
    reader.Update()

    data = reader.GetOutput()
    cells = data.GetPolys()
    triangles = cells.GetData()
    points = data.GetPoints()
    point_data = data.GetPointData()
    Udata = point_data.GetArray(0)

    ntri = triangles.GetNumberOfTuples()/4
    npts = points.GetNumberOfPoints()
    nvls = Udata.GetNumberOfTuples()

    tri = N.zeros((ntri, 3))
    x = N.zeros(npts)
    y = N.zeros(npts)
    ux = N.zeros(nvls)
    uy = N.zeros(nvls)

    for i in xrange(0, ntri):
        tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
        tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
        tri[i, 2] = triangles.GetTuple(4*i + 3)[0]

    for i in xrange(npts):  # modified
        pt = points.GetPoint(i)
        x[i] = pt[0]
        y[i] = pt[1]

    for i in xrange(0, nvls):  # modified
        U = Udata.GetTuple(i)
        ux[i] = U[0]
        uy[i] = U[1]

    return (x, y, tri, ux, uy)

P.clf()
fn = sys.argv[1]  # new
levels = N.linspace(-1.0, 1.0, 21)  # modified
x, y, tri, ux, uy = load_velocity(fn)  # modified
P.tricontour(x, y, tri, ux, levels, linestyles='-',
             colors='black', linewidths=0.5)
P.tricontourf(x, y, tri, ux, levels)  # modified
P.xlim([0, 0.1])  # new
P.ylim([0, 0.1])  # new
P.gca().set_aspect('equal')
P.minorticks_on()
P.gca().set_xticklabels([])
P.gca().set_yticklabels([])
P.title('$U_x$')  # modified
P.savefig("Ux.png ", dpi=300, bbox_inches='tight')
P.savefig("Ux.pdf", bbox_inches='tight')
P.show() # new

Les changements majeurs par rapport au script d'origine sont que le fichier vtk à lire est donné comme argument en utilisant sys.argv, et finalement la figure est affichée à l'écran avec show (). C'est un point.

Exécution du script Python (diagramme de contour de x composante de vitesse)

Le script préparé est destiné à la lecture de données vectorielles, alors exécutons-le sur les données de vitesse de l'échantillon.

$ python plotUx.py postProcessing/surfaces/0.5/U_interpolatedPlane.vtk

Si le diagramme de contour est affiché dans une autre fenêtre sans aucune erreur, il est réussi. Pour quitter, fermez la fenêtre affichée. Puisqu'il s'agit d'un script qui génère également des fichiers png et pdf, vous pouvez confirmer que ces fichiers image sont créés.

Ux.png

Diagramme de contour du vortex

Jusqu'à présent, j'ai dessiné la composante $ x $ de la vitesse, mais je vais également dessiner d'autres quantités physiques.

Calcul du vortex

Au stade où le calcul de la cavité d'OpenFOAM est terminé, la vitesse et la pression sont sorties en tant que grandeurs physiques principales. Dans l'analyse des fluides, le vortex est utilisé pour observer le vortex existant dans la zone d'analyse, alors visualisons également le vortex. Tout d'abord, exécutez l'utilitaire de calcul de vortex standard OpenFOAM vorticity.

$ vorticity -latestTime

Puisque l'option «-latestTime» est ajoutée, le degré de vortex n'est calculé que pour la dernière fois. Les données de «vorticité» sont sorties sous le répertoire «0.5».

Échantillonnage des données de section transversale

Modifiez le sampleDict préparé pour la vitesse et ajoutez les exemples de variables de champ comme suit.

...
fields
(
    p
    U
    vorticity  // new
);
...

Exécutez la commande sample de la même manière que pour la vitesse.

$ sample -latestTime

Après exécution, les données vtk de vorticité (vorticity) seront sorties sous le répertoire postProcessing / surfaces / 0.5.

Préparation du script Python pour vortex

Pour un écoulement de cavité bidimensionnel, le vortex a une composante directionnelle $ z $ perpendiculaire au plan $ x $ - $ y $, il est donc nécessaire de dessiner la composante $ z $ du vecteur vortex. Par conséquent, le script du composant speed $ x $ est étendu comme suit. Nommez le fichier de script «plotWz.py».

plotWz.py


#!/usr/bin/env python

import sys
import os
import numpy as N
import vtk
import matplotlib.pyplot as P

def load_vorticity(filename): # modified
    if not os.path.exists(filename): return None
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.ReadAllVectorsOn()
    reader.Update()

    data = reader.GetOutput()
    cells = data.GetPolys()
    triangles = cells.GetData()
    points = data.GetPoints()
    point_data = data.GetPointData()
    Wdata = point_data.GetArray(0)  # modified

    ntri = triangles.GetNumberOfTuples()/4
    npts = points.GetNumberOfPoints()
    nvls = Wdata.GetNumberOfTuples()  # modified

    tri = N.zeros((ntri, 3))
    x = N.zeros(npts)
    y = N.zeros(npts)
    wx = N.zeros(nvls)  # modified
    wy = N.zeros(nvls)  # modified
    wz = N.zeros(nvls)  # new

    for i in xrange(0, ntri):
        tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
        tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
        tri[i, 2] = triangles.GetTuple(4*i + 3)[0]

    for i in xrange(npts):
        pt = points.GetPoint(i)
        x[i] = pt[0]
        y[i] = pt[1]

    for i in xrange(0, nvls):
        W = Wdata.GetTuple(i)  # modified
        wx[i] = W[0]  # modified
        wy[i] = W[1]  # modified
        wz[i] = W[2]  # new

    return (x, y, tri, wx, wy, wz)  # new

P.clf()
fn = sys.argv[1]
levels = N.linspace(-260, 260, 16)  # modified
x, y, tri, wx, wy, wz = load_vorticity(fn)  # modified
P.tricontour(x, y, tri, wz, levels, linestyles='-',
             colors='black', linewidths=0.5)  # modified
P.tricontourf(x, y, tri, wz, levels)  # modified
P.xlim([0, 0.1])
P.ylim([0, 0.1])
P.gca().set_aspect('equal')
P.minorticks_on()
P.gca().set_xticklabels([])
P.gca().set_yticklabels([])
P.title('$W_z$')  # modified
P.savefig("Wz.png ", dpi=300, bbox_inches='tight')  # modified
P.savefig("Wz.pdf", bbox_inches='tight')  # modified
P.show()

Exécution du script Python (diagramme de contour du composant z du vortex)

Exécutez le script comme dans le cas de la vitesse.

$ python plotWz.py postProcessing/surfaces/0.5/vorticity_interpolatedPlane.vtk

Si le diagramme de contour est affiché dans une autre fenêtre sans aucune erreur, il est réussi.

Diagramme de contour de pression

J'ai dessiné la composante $ x $ de la vitesse et la composante $ z $ du vortex, mais il y a aussi la pression comme quantité physique de sortie. Puisque la pression est une valeur scalaire, il vous suffit de changer la valeur à lire dans le premier script de vitesse de deux à un. La modification du script est simple, je vais donc l'omettre ici. Si vous êtes intéressé, créez-le vous-même. Puisque la pression (p) est déjà spécifiée dans le sampleDict d'OpenFOAM, vous pouvez dessiner un diagramme de contour de la pression en exécutant le script modifié.

Résumé

Dans cet article, j'ai présenté comment visualiser le résultat du calcul d'OpenFOAM sous forme de diagramme de contour avec une coupe bidimensionnelle à l'aide de matplotlib. Selon le système d'exploitation que vous utilisez, il peut être difficile de créer l'environnement Python, mais une fois la construction terminée, vous pouvez modifier les conditions de calcul et générer automatiquement un diagramme de visualisation. Cette fois, l'espacement des lignes de contour (niveaux) est entré manuellement, mais il est également possible de trouver le maximum / minimum et de le déterminer automatiquement sur Python. De plus, matplotlib peut dessiner divers diagrammes de visualisation, veuillez donc le personnaliser vous-même avec les fonctions de Python.

Recommended Posts

Essayez de dessiner un diagramme de contour en utilisant matplotlib avec OpenFOAM
Essayez d'utiliser matplotlib
Dessin graphique avec matplotlib
Essayez d'utiliser matplotlib avec PyCharm
Essayez d'utiliser Spyder inclus dans Anaconda
Essayez d'utiliser LeapMotion avec Python
Essayez d'utiliser l'API Wunderlist en Python
Essayez d'utiliser l'API Kraken avec Python
Essayez d'utiliser la bande HL dans l'ordre
Essayez de dessiner une animation simple en Python
Essayez de dessiner une distribution normale avec matplotlib
Essayez d'utiliser l'API BitFlyer Ligntning en Python
Essayez d'utiliser une classe orientée objet dans R (méthode R6)
Essayez d'utiliser l'API ChatWork et l'API Qiita en Python
Essayez d'utiliser l'API DropBox Core avec Python
Essayez d'utiliser Tkinter
Essayez d'utiliser docker-py
Essayez d'utiliser Cookiecutter
Essayez d'utiliser PDFMiner
Essayez d'utiliser des géopandas
Essayez d'utiliser Selenium
Essayez d'utiliser pandas.DataFrame
Essayez d'utiliser django-swiftbrowser
Essayez d'utiliser tf.metrics
Essayez d'utiliser PyODE
Essayez de dessiner un graphe social à l'aide de l'API Twitter v2
Essayez de supprimer des tweets en masse à l'aide de l'API de Twitter
Essayez la détection des visages en temps réel à l'aide d'une webcam