OpenFOAM演習
pisoFoamによる2次元円柱周りの流れの解析（修正版）
田村
第24回OpenCAE勉強会（岐阜）2013/8/10

1

Reによる流れの変化
双子渦

Re=1.54

wiki.brown.edu
S. Taneda

Re=26

www.dtu.dk
S.Taneda

カルマン渦

Re=1682013/8/10

Re=10000www.aps.org
S.Kumar et al.

www.nptel.iitm.ac.in
Tomas Corke & Hasan Hajib

2

抗力係数

Re=200, CD=1.5
Re=20, CD=2

2013/8/10

流体の力学：中山泰喜著 養賢堂 (1998) P138

3

ストローハル数
http://www.dept.aoe.vt.edu/~jschetz/fluidnature/unit09/unit9a.html

Re=200, St=0.19

The Strouhal number in terms of the Reynolds number of the flow
past a circular cylinder as measured by A. Roshko
2013/8/104NS方程式の数値計算の草分：川口光年さんの計算（1953）

Re=40

週20時間、1年半タイガー計算機で計算。
Kawaguti, M.：J. Phys. Soc. Jpn. 8 (1953), 747.
2013/8/105例題（円柱のCD値を求める）
U=2cm/sp=02次元問題
L1=2cm

L2=4cmy
W=1.5cm

U=2cm/sp=0xO
円柱壁：粘着条件

U=2cm/sp=0

2013/8/106cm

U=2cm/sp=0
W=1.5cm

境界は十分遠方ではないが、今回は時間がないので狭い領域で計算

円柱直径D=1cm 流体は水⇒非圧縮Newton流体
動粘性係数ν=1×10-5m2/s
レイノルズ数Re=UD/ν=20 ⇒層流6支配方程式

OFのデフォルトで条件設定はできない

2013/8/10

7pｉsoFoam
OpenFOAM ver.2.1.x tutorials
incompressible/pisoFoam＜前提＞
・非定常問題・非圧縮性層流・Newton流体/非Newton流体
・非圧縮性乱流モデル（RAS, LES)
＜pisoFoamの由来＞
次の非定常解法を使う
PISO algorithm - Pressure Implicit with Split Operator
2013/8/10

8

演習の流れ円柱stlファイルの作成（FreeCAD）
メッシュ作成Ⅰ（blockMesh）
メッシュ作成Ⅱ（snappyHexMesh）
境界条件、計算条件の設定
計算実行（pisoFoam）
計算結果の表示（paraFoam）
2013/8/10

9円柱のstlファイル作成

2013/8/10

10

円柱作成1（FreeCAD起動）
DEXCSメニューのFreeCADアイコン
をクリックして起動

作業場所：Desktop
2013/8/10

11

円柱作成2（新規作成）

新規作成アイコンをクリック

2013/8/10

起動時の画面12円柱作成3（円柱作成アイコン）

円柱作成アイコンをクリック

2013/8/10

13

円柱作成4（プロジェクトの指定）

プロジェクトをクリック
してハイライト2013/8/10

14

円柱作成5（寸法指定）

半径5mmに変更

Dataシートを選択
2013/8/10

半径5mm、高さ10mmの円柱

15円柱作成6（メッシュ作成）
Meshesメニュ－

デフォルトでOK

プロジェクトを指定

2013/8/10

Meshesメニューの中でCreate mesh from geometryを選択（注：Create mesh from shapeでも可能）
Meshing Toleranceメニューがでる。
メッシュサイズを制御するメニューだがここではデフォルトでOK

16

円柱作成7（stlファイル保存）
ファイル名cylinder.stlを入力

Desktopに保存

Meshを選択
Binary STLで保存（注：ASCIIでも可能）

2013/8/10

Meshesメニューの中でExport Meshを選択
Export meshメニューがでる。
ファイル名 cylinder.stl を指定しBinary STLでDesktopに保存

17

円柱作成8（単位変換）端末アイコン（これを選択）

変換ファイル cylinder-s.stl

元ファイル FreeCADで半径10mmの円柱stlを作ったが単位系はOpenFoamで認識しない
 OpenFoamではmモジュールを使うので1/1000に変換が必要
 端末アイコンをクリック（3つあるが、OpenFOAMコマンドが使えるのは真ん中） Desktopに移動し、以下のコマンドを使って単位変換$cd Desktop$surfaceTransformPoints cylinder.stl cylinder-s.stl -scale ’(0.001 0.001 0.001)’
2013/8/10

18

ベースケースの取得
pisoFoam/ras/cavity

2013/8/10

19

ベースtutorial
ホームフォルダ

コピー

名前の変更 cavity → cylinder

ベースのtutorialのディクショナリをDesktopにコピー
ホーム/OpenFOAM/OpenFOAM-2.1.x/tutorials/incompressibie/pisoFoam/ras/cavity
2013/8/10

ホルダの名称をcavity⇒ cylinder に変更

20ディクショナリ構造の概要
0constantsystem

2013/8/10（時刻ディクショナリ）
U, p など物理量のt=0の時の値、境界条件（constantディクショナリ）
polyMesh/blockMeshDect:領域、メッシュ指定
transportProperties:輸送物性値の指定
turbulenceProperties:乱流モデル（RAS,LES)指定
RASProperties:レイノルズ平均モデルの指定（systemディクショナリ）
controlDict:ソルバー、開始時間、終了時間、タイムステップ等の指定
fvScheme:有限体積法の離散化方法の指定
fvSolution:有限体積法の方程式解法と収束条件の指定
21

0ディクショナリ

U pだけを残す

 0ディクショナリには初期値および境界条件が格納
 U（速度）、 epsilon （乱流エネルギー散逸率）、k（乱流エネルギー）
nuTilda（渦粘性：Spaｌart-Alｌmarasモデルで使用）、 nut（渦粘性）、 p（圧力） U 、p しかここでは使わないので他は削除（置いてあっても問題ない）
 名称を 0 ⇒ 0.org に変更（0のままだとsnappyHexMeshでエラー）
2013/8/1022

BlockMesh作成
基本領域の６面体の設定とメッシュ分割

2013/8/1023

blockMeshDictの格納場所
Cylinderフォルダの中身

blockMeshDict constant/polyMeshの中に基本の6面体メッシュ作成を制御する
blockMeshDictがある。
2013/8/1024

基本領域設定（blockMeshDict）
変更箇所は赤枠

7単位変換 1m⇒0.01m

654
節点の座標 （x y z）

2次元問題とするためz方向厚を1mmと小さくする

23//の後はコメント

0z

1yx
層流ではRe9/4個のセル数が目安
単位長さあたりR3/4！
203/4=9.5なので、円柱直径10mmを10分割する

6面体指定

2013/8/10

x,y,z方向の分割数1セル 1mm角60×30×1=1800 cells
25

基本境界（Patch）設定（blockMeshDict）

7

upANDdown

6
5

4

frontANDback

upstream

3
frontANDback
0z

downstream
2

1

upANDdownyx

2次元問題の時、計算しない面
2013/8/10

26

blockMesh作成

 コマンド実行フォルダへの移動 $cd cylinder
 blockMeshコマンド実行 $blockMesh
2013/8/10

27

blockMesh実行後

変更や追加されたファイル
2013/8/10

28

ParaFoamの起動（BlockMeshの確認）

ハイライトクリック

 ParaFoamの起動
 Applyをクリック
2013/8/10$paraFoam
29

メッシュの確認

Surface With Edgesを指定

2013/8/1030

Patchの確認

Patch名が表示

チェック

2013/8/1031

snappyHexMesh作成
円柱周りのメッシュ作成

2013/8/1032

ベースファイルの取得

 tutorials/incompressible/pisoFoam/les/motorBike/motorBike/systemから以下の3ファイルをコピー①snappyHexMeshDict （snappyHexMeshの制御ファイル）②decomposeDict （並列計算で使用。これがないとsnappyでエラー）③forceCoeffs（Cd値計算で使用。ついでにコピー）
2013/8/10

33

snappyHexDictの変更
１段階（○）：表面に適合したメッシュ除去
２段階（○）：表面メッシュのスナップ（平滑化）
３段階（×）；表面へのレイヤー挿入

円柱ファイルの名称
ファイルタイプ：stl
表面の名前指定：cylinder

2013/8/1034

snappyHexMeshDictの変更
細分化したい表面の名称
細分化の程度 （最小、最大）数が0⇒1⇒2⇒・・・大きくなるほど細分化
最小レベルは全ての表面に適用
最大レベルはresolveFeatureAngel 以上の角度の表面に適用細分化したいBox領域の細分化
今回使わないのでコメントアウト

Snappyを適用する領域を指定。今回は円柱の外部の座標を選ぶ。セルの表面の座標は好ましくない。
但し、今回の点はセル表面であるが、認識できた。

2013/8/10

35

snappyHexMesh実行 snappyHexMesh実行2013/8/10$snappyHexMesh
36

snappyHexMesh実行 snappy成功 cell数は1720
 0.005、0.01ディレクトリが生成（後に説明するcotrolDictの時間ディレクトリ作成間隔が0.005secの場合）

2013/8/10

37

メッシュの確認（第1段階）
次の時間に進める

時間 0.005

 $paraFoamを起動
 円柱部分のセルが除去されている
2013/8/10
 表面は階段状

38

メッシュの確認（第2段階）
時間 0.01

 円柱表面の平滑化
2013/8/10

39

メッシュ確認

 円柱側面は1層だけ切れている。（2次元問題OK）
 円柱側面のPatchはcylinder_zone0の名前が付けられている。
2013/8/10

40

snappyHexMesh(再)0.005 0.01 を削除

 生成した0.01ディレクトリのメッシュを使って0ディレクトリを作っても良いが別の方法をトライ。一回で0ディレクトリが生成 $snappyHexMesh –overwrite
2013/8/10
 $rm -r 0.005 0.01 (0.005 0.01を削除）

41

0ディクショナリ作成
初期値・境界値の設定

2013/8/10

42

U
constant/polyMesh/boundary（Patch名とtype等を記述）

0.org/U次元(kg m s K mol A cd)=(m/s)
領域初期値 全域2cm/s

固定値 上流2cm/s
チュートリアルでよく使っているで使用流出条件（流入あり）
流入時 速度2cm/s
流出時 内部メッシュ値
流出条件（流入あり）
流入時 速度2cm/s
流出時 内部メッシュ値

2次元問題条件

2013/8/10

固定値 粘着条件 43

p=(p/ρ)次元=(m2/s2)0.org/p
初期値 全域 0

p=0

p=0

p=0

2次元問題条件

2013/8/10

あまり適切ではないが、選択できる条件では、これが最もよい（付録参照）

440ディクショナリ作成

コピー 0ディクショナリにはU、pファイルがないので0.orgからコピー
2013/8/10

45

Constantディクショナリのその他ファイルの設定
物性値と流体モデル選定

2013/8/1046物性値と解析モデル
turbulencePropertiesファイル（変更なし）

RASPropertiesファイル

乱流モデルの選択
RASmodel （Reynolds平均モデル）
LESmodel （Large Eddy Simulation モデル）
tranportPropertiesファイル（変更なし）
動粘性係数[m2/s]

輸送係数モデル選択
Newtonian （ニュートン流体）
CrossPowerLawCoefss （非ニュートン流体）など
2013/8/10

層流47

systemディクショナリのファイル設定
計算スキームと計算制御の設定

2013/8/1048

controlDict設定⇒計算ソルバー⇒startTimeの設定時間から計算開始⇒計算開始時間⇒endTimeの設定時間に計算終了⇒計算終了時間（試行錯誤で決定）⇒計算時間ステップ（クーラン数で決定）⇒witeIntervalごとに時刻ディクショナリ生成⇒1000ステップ、1秒ごと⇒時刻ディレクトリ最大数：制限なし⇒データファイルフォーマット選択⇒writeFormatでの有効桁数⇒データ圧縮指定⇒timeFormatのフォーマット選択⇒timeFormat の設定で用いる指数⇒各ディクショナリのyes/noスイッチ⇒物体に係る力の係数計算用の関数
2013/8/1049

forceCoeffs

円柱Patch密度[kg/m3]
揚力方向 y抗力方向 x
ピッチモーメントの回転軸 z
上流での流速 [m/s]

物体中心
代表寸法 [m]射影面積 [m2]2013/8/1050

fvscheme (変更なし)

2013/8/1051

ｆｖscheme (変更なし)

2013/8/1052

Fvsolution(変更なし)

2013/8/1053

STLファイルの格納

2013/8/1054

triSurfaceへ格納

 snappyHexMeshで指定したcylinder-s.stlについて
constant/triSurfaceディクショナリを見に行く constantにtriSurfaceを作成
2013/8/10
 Desktopのcylinder-s.stlをtriSurfaceに移動

55

計算実行

2013/8/10

56

pisoFoamの実行 $pyFoamPlotRunner.py pisoFoam
 pyFoamPlotRunner.pyで残差等がモニターできる。
2013/8/10

57

計算終了時の画面

計算時間
抗力係数 （実験値は2.0程度）
揚力係数
ピッチモーメント係数

2013/8/1058

Cdの履歴10秒程度で収束
2013/8/1059

paraFoam
スケールをクリック

2013/8/10

Uを選択

時間経過を観察

図の拡大⇒マウス中軸を回す

60

流線作図
StreamTracer クリック

クリック

2013/8/10

61

流線作図

⑤クリック

①Line Sourceを選択③0.007を入力②Y Axisを選択④50を入力2013/8/10

62

流線作図

2013/8/1063速度分布と流線の合図

目玉アイコンをON

Opacity（透明度）
1.00⇒0.40

2013/8/10

双子渦もそれらしい64

まとめ
pisoFoamを用いて2次元円柱周り流れの解析の演習を実施
解析の流れは以下のとおり。
①円柱STLファイル作成（FreeCAD)②blockMeshファイル作成・実行③snappyHexMeshファイル作成・実行④各種ファイル(境界、物性、制御）の作成⑤pisoFoamの実行⑥paraFoamを使った流線図の作成
2013/8/1065

付録Ⅰ
境界条件、メッシュ影響

2013/8/1066

境界条件pの影響
PatchUp(case1)p(case2)p(case3)p(case4)Mesh-

60×30×1←←←

upstream(2,0,0)00

downstream(2,0,0)00

upANDdown(2,0,0)00cylinder_zone0(0,0,0)zeroGradient

fixedValue;value$internalField

fixedValue;value$internalField

zeroGradientfrontANDback

empty

←

←

←

←

Cd-

4.13

3.61

4.35

5.00

双子渦-

△××△

計算時間-

136

136

139

150

2013/8/10

zeroGradient zeroGradient00

zeroGradient zeroGradient Cd値は実験（Cd=2.0)とかなり異なる 特に、円筒壁には＄internalFieldを使うのは不適

67境界条件pの影響

2013/8/10

case1

case2

case3

case4

双子渦の様子が実験と異なる

68

計算領域の影響
PatchUp(case1)p(case5)p(case6)

Δ x-Δ y-6-3cm

24-12cm

28-16cmBase Mesh-

60×30

240×120280×160

upstream(2,0,0)0←←

downstream(2,0,0)0←←upANDdown(2,0,0)0←←cylinder_zone0(0,0,0)zeroGradient←←frontANDback

empty

empty

←

←

Cd-

4.13

2.26

2.20

双子渦-

△○○

計算時間

-

136

1,404

2,054

2013/8/10

 4倍程度、領域が広くなるとCdが収束してくる 但し計算時間は10倍69計算領域の影響(case5)

2013/8/1070

付録Ⅱ
Re=200の場合

2013/8/1071

領域とメッシュの影響
項目
Δ x-Δ yBase MeshCells
CdSt
計算時間(sec)

case1
6-3cm
60×30
1.06
0.195
129(2.2min)

case2
case3
24-12cm 24-12cm
240×120 480×240
114,884
1.19
1.36
0.186
0.186
2,729
19,987(45.5min) (5.5hrs)

実験ではCd=1.5、St=0.19なのでCase3では程々の精度
2013/8/10

72

CdとClの履歴(case3）

2013/8/10

73

流れの様子（case3）

2013/8/10

74

付録Ⅲ
extrudeMeshを使った
2Dメッシュの作成

2013/8/10

75

概要
オープンCAE勉強会＠富山/20130525第10回での富山大学の中川先生の資料を参考
http://eddy.pu-toyama.ac.jp/オープンCAE勉強会-富山/#_99

snappyHexMeshで3次元メッシュを作成。extrudeMeshでPatch表面メッシュだけを押し出して2次元メッシュとする。
演習での方法では円柱周りメッシュが１層にできない場合があるが、本方法では必ず１層にできる。
本方法で領域のメッシュを部分的に詳細にすることが可能で計算時間の短縮になる。
2013/8/1076blockMeshDictの修正

分解

2013/8/10

77

snappyHexMeshDictの修正1＜詳細メッシュ領域の指定＞
２つのボックスを指定
refinementBox1
refinementBox2（名前は任意）
領域名称

指定図形タイプ：ボックスボックスを規定できる2点を指定

領域名称

指定図形タイプ：ボックスボックスを規定できる2点を指定
2013/8/10

78

snappyHexMeshDictの修正2＜詳細メッシュ領域の分割法指定＞
領域名称

ボックス内部／外部を指定
細分化レベル levels（距離、レベル）
領域名称

ボックス内部／外部を指定
細分化レベル levels（距離、レベル）
refinementBox1の細分化レベルが1に対しrefinementBox2はレベル2でより細かく細分化される
2013/8/10

79extrudeMeshDictの作成
pimpleDyMFoam/wingMotion/wingMotion2D_simpleF
oam/system/extrudeMeshDictを/systemにコピー
現在のディレクトリでソースケースを参照

押し出すPatch： front

押し出されたPatch： back押し出す層数2013/8/10

80

Uとpの修正

frontANDbackの分離

2013/8/10

81

snappyHexMesh後のメッシュ$blockMesh$snappyHexMesh -overwrite

細かい領域だけ5層

2013/8/10

82

extrudeMeshの実行後のメッシュ$extrudeMesh

2013/8/10

一層だけ

83

CdとStの比較
項目
Δ x-Δ yCells
CdSt
計算時間(sec)

2013/8/10

今回
24-12cm
38,5561.34
0.195
1,880(31.3min)

等間隔メッシュ
24-12cm
114,8841.360.186
19,987(5.5hrs)84

流線
等間隔メッシュ

今回メッシュ

2013/8/1085