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

18568 ワード
Python matplotlib vtk OpenFOAM OpenCAE Python テキストリンク

はじめに

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を使って圧力と流速コンター図を描くとこんな感じになります。

データのサンプリングには 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つ目の引数を調整すると良いです。

流線が必要な場合は 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 で矢印の形状を指定できます。

まとめ

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

以上です。