Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

OpenFOAMの結果をVTKで出力してmatplotlibで表示する

はじめに

OpenFOAMから得られる圧力や流速などのデータを,paraViewを介さずにmatplotlibを使って2次元のコンター図として表示する方法についてまとめます。

解析の度にparaFoamを使ってコンター図に流線を追加したり,カラーバーを調整するのは時間がかかってしまいますので,postProcessコマンドのsurfaces を使って計算結果をVTK形式で出力させ,線上や面上のデータを取得してmatplotlibで表示させてみます。

実験した環境は以下の通りです。

  • Ubuntu 18.04LTS
  • OpenFOAM-v1912
  • Python 3.6.9
  • matplotlib 3.3.3

準備

VTKをPythonで扱うため,pip3を使ってインストールしておきます。

$ pip3 install vtk

データの用意

データとして,icoFoam のチュートリアルケース cavity を使います。

$ cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
$ cd cavity
$ blockMesh
$ icoFoam

paraFoamを使って圧力と流速コンター図を描くとこんな感じになります。
contour.png

データのサンプリングには postProcess コマンドの surfaces を使います。そのために system/surfaces を用意します。

system/surfaces
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
-------------------------------------------------------------------------------
Description
    Writes out surface files with interpolated field data in VTK format, e.g.
    cutting planes, iso-surfaces and patch boundary surfaces.

    This file includes a selection of example surfaces, each of which the user
    should configure and/or remove.

\*---------------------------------------------------------------------------*/

#includeEtc "caseDicts/postProcessing/visualization/surfaces.cfg"

fields       (U p);
surfaceFormat  vtk;

surfaces
(
    xy  
    {
        type                 plane;
        planeType            pointAndNormal;
        pointAndNormalDict
        {
            point       (0 0 0.005);
            normal    (0 0 1);
        }
    }
);

// ************************************************************************* //

fields で速度と圧力のデータを,surfaceFormat で vtk を指定しています。
平面のデータを取得するために,surfaces の type で plane を指定します。pointAndNormal は基準点と法線ベクトルの向きを指定します。
cavity のデータは z 軸方向に0.01の厚みがありますので,ここでは厚みの中心を通る xy 平面を抽出するように z=0.005 としています。

postProcess の実行

surfaces オプションをつけて実行します。

$ postProcess -func surfaces

$FOAM_RUN/cavity/postProcessing/surfaces の各時刻毎のディレクトリ内に vtp ファイルが作られます。
vtk と指定したはずなのに vtp 形式で出力されて焦りますが,落ち着いて vtk_to_numpy を使って読み込んでいきます。

matplotlib でコンター図出力

vtp 形式の polygon データを読み込む場合は vtkXMLPolyDataReader() を使います。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from scipy.interpolate import griddata


filename = "postProcessing/surfaces/0.5/xy.vtp"

reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(filename)
reader.Update()
data = reader.GetOutput()

# cell data から point data変換
cell2point = vtk.vtkCellDataToPointData()
cell2point.SetInputData(reader.GetOutput())
cell2point.Update()

# 座標データの配列化
points = data.GetPoints()
coord = vtk_to_numpy(points.GetData()) # (x,y,z)座標の2次元配列
x = coord[:,0]
y = coord[:,1]
z = coord[:,2] 

# メッシュグリッド用
xmin, xmax = min(x), max(x)
ymin, ymax = min(y), max(y)
xi = np.linspace(xmin, xmax, 100)
yi = np.linspace(ymin, ymax, 100)

# GetAbstractArray(0)は圧力、GetAbstractArray(1)は速度ベクトルデータ
p = vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray(0))
u = vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray(1))
speed = np.sqrt(u[:,0]**2 + u[:,1]**2) # ベクトルをスカラー値に変換

# 圧力のコンター図出力
levels = np.linspace(-5,5,21) # 描画する範囲の設定(最小,最大,色の分割数)
plt.tricontourf(x,y,p,levels=levels,cmap="jet",vmin=-5, vmax=5) #vmin,vmax で色の範囲を設定
cbar = plt.colorbar()
cbar.set_ticks(np.arange(-5,5.1,0.5)) # カラーバーの目盛り
cbar.set_label("p")
plt.show()

# 流速のコンター図出力
levels = np.linspace(0,1,21) 
plt.tricontourf(x,y,speed,levels=levels,cmap="jet", vmin=0, vmax=1)
cbar = plt.colorbar()
cbar.set_ticks(np.arange(0,1.01,0.1)) 
cbar.set_label("U")
plt.show()

このスクリプトを実行して出力される圧力と流速グラフがこちら。
色の変化をもっと滑らかにしたい場合は,linspace の3つ目の引数を調整すると良いです。

pU.png

流線が必要な場合は streamplot() を使って描画することができます。

# 流線出力
levels = np.linspace(0,1,21) 
velocity = griddata((x, y), u, (xi[None,:], yi[:,None]), method='linear') # 速度ベクトルをグリッドに合わせて線形補間
speed2 = np.sqrt(velocity[:,:,0]**2 + velocity[:,:,1]**2)
strm = plt.streamplot(xi, yi, velocity[:,:,0], velocity[:,:,1], linewidth=1, arrowstyle="-", density=1.5)
plt.tricontourf(x,y,speed,levels=levels,cmap="jet", vmin=0, vmax=1)
cbar = plt.colorbar()
cbar.set_ticks(np.arange(0,1.01,0.1)) 
cbar.set_label("U")
plt.show()

density で流線の密度,arrowstyle で矢印の形状を指定できます。
stream.png

まとめ

同じケースで何度も計算をする場合,いちいちparaFoamを起動して操作する手間が省けるので便利です。
サクッとコンター図だけでも確認したい場合に重宝します。

以上です。

hiro_kuramoto
オープンソースソフトを使った多目的最適化に取り組んでいます。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away