[PYTHON] 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.
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