3
Help us understand the problem. What are the problem?

posted at

【スマホでCFD】OpenFOAM@Google Colaboratoryで可視化まで行う

はじめに

Google ColabでOpenFOAMチュートリアルを実行し、計算結果の可視化まで行う方法を記載する。

以下の記事でColabでOpenFOAMを実行する方法が紹介されているが、本記事では可視化までColabで行う。

【最速】5分で実施OpenFOAM@Google Colaboratory

OpenFOAM@Colabの利点

  • クラウドの壊しても問題ない環境で環境構築できる
  • PCが無い場合でも、スマホでOpenFOAMを体験できる
    スマホのブラウザでColabノートブックにアクセスすると、スマホで計算を実行し、その場で可視化まで可能

OpenFOAM@Colabの実行

Colabにアクセスし、ノートブックを新規作成。

1. OpenFOAMのインストール

下記をセルに入力し、ColabのサーバにOpenFOAMをインストール。

%%bash
sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -"
add-apt-repository http://dl.openfoam.org/ubuntu
apt-get update
apt-get -y install openfoam9
pip3 install vtk #可視化のためにvtkをインストール

2. 計算Caseを準備

チュートリアルケースをカレントディレクトリにコピーする。
ここではCavityケースを利用する。

%%bash
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity . #チュートリアルケースをコピー
%cd ./cavity

3. 計算実行

下記セルを実行し、計算を実行する。

%%bash
. /opt/openfoam9/etc/bashrc #環境変数の読み込み
blockMesh >log.blockMesh 2>&1 #メッシュ作成
icoFoam >log.icoFoam 2>&1 #ソルバー実行
ls #計算結果が出力されていることを確認

CavityResult.png

4. 可視化

以下セルで、計算結果をvtkファイルに変換する。

%%bash
. /opt/openfoam9/etc/bashrc
foamToVTK #計算結果をVTK形式に変換

以下セルで、こちらを参考に、matplotlibでvtkファイルを可視化する関数を定義する。

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
from ipywidgets import interact

def plotVtk(t):
  filename = f"VTK/cavity_{t}.vtk"

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

  # cell data から point data変換
  cell2point = vtk.vtkCellDataToPointData()
  cell2point.SetInputData(reader.GetOutput())
  cell2point.Update()
  # 座標データの配列化
  coord = vtk_to_numpy(cell2point.GetOutput().GetPoints().GetData())

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

  # 圧力のコンター図出力
  plt.axes().set_aspect('equal')
  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()

  # 流速のコンター図出力
  plt.axes().set_aspect('equal')
  levels = np.linspace(-1,1,21) 
  plt.tricontourf(x,y,U,levels=levels,cmap="jet", vmin=-1, vmax=1)
  cbar = plt.colorbar()
  cbar.set_ticks(np.arange(-1.0,1.01,0.1)) 
  cbar.set_label("u")
  # plt.quiver(x,y,U,V,color='red',angles='xy',scale_units='xy', scale=10.0) #長さでベクトル場を可視化
  plt.quiver(x,y,U/speed,V/speed,speed,cmap='jet',scale=15.0) #カラーマップでベクトル場を可視化
  plt.show()

  # 流線出力
  plt.axes().set_aspect('equal')
  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()

t=100時点を可視化する。

plotVtk(t=100)

Cavity100p.png
Cavity100u.png
Cavity100stream.png

ノートブックなので、インタラクティブな表示もできる。

interact(plotVtk, t=(0,100,20) )

CavityInteractive.png

実際のノートブックはこちら

拡張性

簡単な計算条件の変更はsedコマンド等で可能。

!sed -i -e "25c \        value           uniform (2 0 0);" 0/U #上壁面の境界条件を2m/sに変更

まとめ

Google ColabでOpenFOAMを実行し、matplotlibで可視化を行った。
Google ColabではGUIをサポートしていないらしく、直接PawaviewやparaFoamが利用できないため、かなり面倒な方法になった。
Paraviewは偉大。

おまけ

ダムブレークのチュートリアルも同様にOpenFOAM@Colabで可視化することができた。
DamBreakのノートブック
damBreak.gif

参考

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
Sign upLogin
3
Help us understand the problem. What are the problem?