소개

OpenFOAM에서 얻은 압력과 유속 등의 데이터를 paraView를 통하지 않고 matplotlib을 사용하여 2 차원의 등고선으로 표시하는 방법을 요약합니다.

해석마다 paraFoam을 사용하여 등고선 다이어그램에 유선을 추가하거나 색상 막대를 조정하는 데 시간이 걸리므로 postProcess 명령의 surfaces를 사용하여 계산 결과를 VTK 형식으로 출력하여 선상 혹은 면상의 데이터를 취득해 matplotlib로 표시시켜 보겠습니다.

실험한 환경은 다음과 같습니다.

우분투 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의 세 번째 인수를 조정하는 것이 좋습니다.



유선이 필요한 경우는 streamplot() 를 사용해 paint 할 수 있습니다.

# 流線出力
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을 기동해 조작하는 수고를 줄일 수 있으므로 편리합니다.
사쿠와 등고선 그림만으로도 확인하고 싶은 경우에 편리합니다.

이상입니다.