[PYTHON] Versuchen Sie, ein Konturdiagramm mit matplotlib mit OpenFOAM zu zeichnen

Einführung

In diesem Artikel geht es um Tipps von OpenFOAM, einer Open-Source-Software für die Flüssigkeitsanalyse. In OpenFOAM wird das Berechnungsergebnis mit ParaView als Nachbearbeitung (Visualisierung) gezeichnet, das gezeichnete Bild ist jedoch grundsätzlich im Rasterformat und zweidimensional. Die Kontur des Querschnitts kann nicht im Vektorformat ausgegeben werden. Wenn der zu visualisierende Querschnitt / die zu visualisierende Variable festgelegt wird, wird das Starten von ParaView und das Zeichnen mühsam. In diesem Artikel wird daher zusammengefasst, wie eine zweidimensionale Schnittkontur halbautomatisch gezeichnet und mit matplotlib, einer Visualisierungsbibliothek von Python, im Vektorformat ausgegeben wird.

Umgebung

Die diesmal erforderliche Software ist OpenFOAM und Python. Die Version von OpenFOAM wurde kürzlich auf 4.0 oder v1606 aktualisiert. Dieser Artikel konzentriert sich jedoch auf 2.4.x. Für Python ist 2.7 das Ziel. Numpy, matplotlib, vtk (python-vtk), tkinter (python-tk) sind als Python-Pakete erforderlich. Wenn nicht, installieren Sie sie bitte zusätzlich.

Überprüfen Sie nach dem Erstellen der Umgebung den Betrieb mit dem Terminal.

OpenFOAM

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

foamInstallationTest überprüft die grundlegenden Umgebungseinstellungen des installierten OpenFOAM. In den folgenden Befehlen wird der Hohlraumfluss ausgeführt, der das grundlegende Tutorial von OpenFOAM darstellt. Wenn bisher keine Fehler vorliegen, ist dies in Ordnung.

Python

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

Wenn jedes Paket (Modul) installiert ist, wird nichts angezeigt und die Eingabeaufforderung wird normal zurückgegeben. Wenn der Fehler "ImportError: Kein Modul mit dem Namen xxxx" zurückgegeben wird, fügen Sie das entsprechende Paket erneut ein.

Erfassung zweidimensionaler Querschnittsdaten (vtk)

Erfassen Sie zunächst aus dem Analyseergebnis von OpenFOAM die zweidimensionalen Querschnittsdaten, die Sie im vtk-Format zeichnen möchten. Sie können OpenFOAM-Daten mit ParaView lesen und Querschnittsdaten mit Slice ausgeben. In dieser Zeit werden Konturdiagramme jedoch so oft wie möglich automatisch mit Befehlen oder Skripten gezeichnet, sodass ParaView eine GUI-Operation erfordert Nicht benutzt. Hier verwenden wir das OpenFOAM-Dienstprogramm "sample" (das in "postProcess" in OpenFOAM 4.0 und höher integriert ist).

Erstellen Sie Querschnittsdaten mit den Ergebnissen des Hohlraum-Tutorials, das bei der Überprüfung der Umgebungseinstellungen durchgeführt wurde. Um sample zu verwenden, kopieren und bearbeiten Sie zunächst die Vorlage von sampleDict, der Einstellungsdatei dieses Befehls.

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

Die ursprüngliche Vorlagendatei beschreibt verschiedene Stichprobenverfahren, einschließlich Verwendungskommentare. Dieses Mal werden wir jedoch nur die $ x $ - $ y $ -Ebene des zweidimensionalen Hohlraums abtasten, also einfach "sampleDict" wie folgt. Zu

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

Da sets beim Erfassen eindimensionaler Daten wie z. B. geraden Linien verwendet wird, lassen Sie es diesmal leer. Es gibt verschiedene Möglichkeiten, "Oberflächen" anzugeben, aber dieses Mal verwenden wir die durch Interpolation erhaltene.

Führen Sie nach dem Bearbeiten von "sampleDict" den Befehl "sample" aus, um zweidimensionale Querschnittsdaten zu generieren. Mit der Option "-latestTime" werden nur die Daten des letzten Males (0,5) abgetastet.

$ sample -latestTime

Wenn der Befehl ausgeführt wird, wird im aktuellen Verzeichnis ein Verzeichnis mit dem Namen "postProcessing / surface" erstellt und zum angegebenen Zeitpunkt vtk-Daten ausgegeben.

Python-Skriptvorbereitung

Das Format der vtk-Datei, die mit dem oben genannten Dienstprogramm "sample" erstellt wurde, ist das sogenannte Legacy-Format Version 2.0. Wenn Sie in Python Daten im vtk-Format lesen und an matplotlib der Visualisierungsbibliothek übergeben können, können Sie frei zeichnen. Sie können das Skript selbst schreiben, aber Alexey Matveichev hat [Plotten von VTK-Dateien mit matplotlib und python-vtk] geschrieben (http://matveichev.blogspot.jp/2014/03/plotting-vtk-files-with-matplotlib] -and.html) ist geschrieben, verwenden Sie also das veröffentlichte Skript. Das im Artikel veröffentlichte Skript enthält eine vtk-Version von python-vtk der 5.x-Serie. Wenn Sie also eine Umgebung mit der neuen 6.x-Serie erstellen, wird es im Artikel [auf github veröffentlichtes Skript](https :: Bitte verwenden Sie //gist.github.com/mrklein/44a392a01fa3e0181972). Hier wird davon ausgegangen, dass das 5.x-System verwendet wird.

Kopieren Sie das ursprüngliche Skript, legen Sie es im OpenFOAM-Parsing-Verzeichnis ab (nennen wir es "plotUx.py") und ändern Sie es ein wenig.

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

Die wichtigsten Änderungen gegenüber dem ursprünglichen Skript bestehen darin, dass die zu lesende vtk-Datei mit sys.argv als Argument angegeben wird und schließlich die Abbildung mit show () auf dem Bildschirm angezeigt wird. Es ist ein Punkt.

Ausführung des Python-Skripts (Konturdiagramm der x-Geschwindigkeitskomponente)

Das vorbereitete Skript dient zum Lesen von Vektordaten. Führen Sie es also für die Daten der Abtastgeschwindigkeit aus.

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

Wenn das Konturdiagramm fehlerfrei in einem anderen Fenster angezeigt wird, ist es erfolgreich. Schließen Sie zum Beenden das angezeigte Fenster. Da es sich um ein Skript handelt, das auch PNG und PDF ausgibt, können Sie bestätigen, dass diese Bilddateien erstellt wurden.

Konturdiagramm des Wirbels

Bisher habe ich die $ x $ -Komponente der Geschwindigkeit gezeichnet, aber ich werde auch andere physikalische Größen zeichnen.

Berechnung des Wirbels

In dem Stadium, in dem die Hohlraumberechnung von OpenFOAM abgeschlossen ist, werden Geschwindigkeit und Druck als physikalische Hauptgrößen ausgegeben. Bei der Flüssigkeitsanalyse wird der Wirbel verwendet, um die im Analysebereich vorhandenen Wirbel zu beobachten. Visualisieren wir also auch den Wirbel. Führen Sie zunächst das OpenFOAM-Standarddienstprogramm zur Wirbelberechnung "Vorticity" aus.

$ vorticity -latestTime

Da die Option "-latestTime" hinzugefügt wird, wird der Wirbel nur zum letzten Mal berechnet. Die "Vorticity" -Daten werden im "0.5" -Verzeichnis ausgegeben.

Abtastung von Querschnittsdaten

Bearbeiten Sie das für die Geschwindigkeit vorbereitete "sampleDict" und fügen Sie die Beispielfeldvariablen wie folgt hinzu.

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

Führen Sie den Befehl sample wie bei der Geschwindigkeit aus.

$ sample -latestTime

Nach der Ausführung werden vtk-Vorticity-Daten ("Vorticity") im Verzeichnis "postProcessing / surface / 0.5" ausgegeben.

Vorbereitung des Python-Skripts für den Wirbel

Für eine zweidimensionale Hohlraumströmung hat der Wirbel eine $ z $ -Richtungskomponente senkrecht zur $ x $ - $ y $ -Ebene, daher ist es notwendig, die $ z $ -Komponente des Wirbelvektors zu zeichnen. Daher wird das Skript für die Komponente speed $ x $ wie folgt erweitert. Nennen Sie die Skriptdatei "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()

Ausführung des Python-Skripts (Konturdiagramm der z-Komponente des Wirbelgrades)

Führen Sie das Skript wie im Fall der Geschwindigkeit aus.

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

Wenn das Konturdiagramm fehlerfrei in einem anderen Fenster angezeigt wird, ist es erfolgreich.

Druckkonturdiagramm

Ich habe die $ x $ -Komponente der Geschwindigkeit und die $ z $ -Komponente des Wirbels gezeichnet, aber es gibt auch Druck als physikalische Ausgangsgröße. Da der Druck ein skalarer Wert ist, müssen Sie nur den im ersten Geschwindigkeitsskript zu lesenden Wert von zwei auf eins ändern. Das Ändern des Skripts ist einfach, daher werde ich es hier weglassen. Wenn Sie interessiert sind, erstellen Sie es bitte selbst. Da der Druck (p) bereits im sampleDict von OpenFOAM angegeben ist, können Sie ein Konturdiagramm des Drucks zeichnen, indem Sie das geänderte Skript ausführen.

Zusammenfassung

In diesem Artikel habe ich vorgestellt, wie das Berechnungsergebnis von OpenFOAM als Konturdiagramm mit zweidimensionalem Querschnitt mithilfe von matplotlib visualisiert wird. Abhängig vom verwendeten Betriebssystem kann es schwierig sein, die Python-Umgebung zu erstellen. Sobald die Erstellung abgeschlossen ist, können Sie die Berechnungsbedingungen ändern und automatisch ein Visualisierungsdiagramm erstellen. Dieses Mal wird der Konturlinienabstand (Ebenen) manuell eingegeben, es ist jedoch auch möglich, das Maximum / Minimum zu finden und es automatisch in Python zu bestimmen. Darüber hinaus kann matplotlib verschiedene Visualisierungsdiagramme zeichnen. Passen Sie diese daher zusammen mit den Funktionen von Python selbst an.