PisoFoamによる　　　　　　　　　2次元円柱周りの流れの解析

Size: px

Start display at page:

Download "PisoFoamによる　　　　　　　　　2次元円柱周りの流れの解析"

なおよどぎみ
5 years ago
Views:

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

2 Re による流れの変化双子渦 Re=1.54 wiki.brown.edu S. Taneda Re=26 S.Taneda カルマン渦 Re=168 Re= /8/ S.Kumar et al. Tomas Corke & Hasan Hajib

3 抗力係数 Re=200, C D =1.5 Re=20, C D =2 流体の力学 : 中山泰喜著養賢堂 (1998) P /8/10 3

4 ストローハル数 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/10 4

5 NS 方程式の数値計算の草分 : 川口光年さんの計算 (1953) Re=40 週 20 時間 1 年半タイガー計算機で計算 Kawaguti, M.:J. Phys. Soc. Jpn. 8 (1953), /8/10 5

6 例題 ( 円柱の C D 値を求める ) 2 次元問題 L 1 =2cm U=2cm/s p=0 L 2 =4cm y W=1.5cm U=2cm/s p=0 O x U=2cm/s p=0 円柱壁 : 粘着条件 W=1.5cm U=2cm/s p=0 6cm 境界は十分遠方ではないが今回は時間がないので狭い領域で計算円柱直径 D=1cm 流体は水非圧縮 Newton 流体動粘性係数 ν= m 2 /s レイノルズ数 Re=UD/ν=20 層流 2013/8/10 6

7 支配方程式 OF のデフォルトで条件設定はできない 2013/8/10 7

8 pisofoam OpenFOAM ver.2.1.x tutorials incompressible/pisofoam < 前提 > 非定常問題非圧縮性層流 Newton 流体 / 非 Newton 流体非圧縮性乱流モデル (RAS, LES) <pisofoam の由来 > 次の非定常解法を使う PISO algorithm - Pressure Implicit with Split Operator 2013/8/10 8

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

10 円柱の stl ファイル作成 2013/8/10 10

11 円柱作成 1(FreeCAD 起動 ) DEXCS メニューの FreeCAD アイコンをクリックして起動作業場所 :Desktop 2013/8/10 11

12 円柱作成 2( 新規作成 ) 新規作成アイコンをクリック起動時の画面 2013/8/10 12

13 円柱作成 3( 円柱作成アイコン ) 円柱作成アイコンをクリック 2013/8/10 13

14 円柱作成 4( プロジェクトの指定 ) プロジェクトをクリックしてハイライト 2013/8/10 14

15 円柱作成 5( 寸法指定 ) 半径 5mmに変更 Dataシートを選択半径 5mm 高さ 10mm の円柱 2013/8/10 15

16 円柱作成 6( メッシュ作成 ) Meshes メニュ - プロジェクトを指定デフォルトで OK Meshesメニューの中でCreate mesh from geometryを選択 ( 注 :Create mesh from shapeでも可能 ) Meshing Toleranceメニューがでるメッシュサイズを制御するメニューだがここではデフォルトでOK 2013/8/10 16

17 円柱作成 7(stl ファイル保存 ) ファイル名 cylinder.stl を入力 Desktop に保存 Mesh を選択 Binary STL で保存 ( 注 :ASCII でも可能 ) Meshesメニューの中でExport Meshを選択 Export meshメニューがでるファイル名 cylinder.stl を指定しBinary STLでDesktopに保存 2013/8/10 17

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

19 ベースケースの取得 pisofoam/ras/cavity 2013/8/10 19

20 ベース tutorial ホームフォルダコピー名前の変更 cavity cylinder ベースの tutorial のディクショナリを Desktop にコピーホーム /OpenFOAM/OpenFOAM-2.1.x/tutorials /incompressibie/pisofoam/ras/cavity ホルダの名称を cavity cylinder に変更 2013/8/10 20

ディクショナリ構造の概要 0 constant system ( 時刻ディクショナリ ) U, p など物理量の t=0 の時の値境界条件 (constant ディクショナリ ) polymesh/blockmeshdect: 領域メッシュ指定 transportproperties: 輸送物性値の指定

21 ディクショナリ構造の概要 0 constant system ( 時刻ディクショナリ ) U, p など物理量の t=0 の時の値境界条件 (constant ディクショナリ ) polymesh/blockmeshdect: 領域メッシュ指定 transportproperties: 輸送物性値の指定 turbulenceproperties: 乱流モデル (RAS,LES) 指定 RASProperties: レイノルズ平均モデルの指定 (system ディクショナリ ) controldict: ソルバー開始時間終了時間タイムステップ等の指定 fvscheme: 有限体積法の離散化方法の指定 fvsolution: 有限体積法の方程式解法と収束条件の指定 2013/8/10 21

モデルで使用 ) nut( 渦粘性 ) p ( 圧力 ) U p しかここでは使わないので他は削除 (

22 0 ディクショナリ U p だけを残す 0 ディクショナリには初期値および境界条件が格納 U( 速度 ) epsilon ( 乱流エネルギー散逸率 ) k( 乱流エネルギー ) nutilda( 渦粘性 :Spalart-Allmaras モデルで使用 ) nut( 渦粘性 ) p ( 圧力 ) U p しかここでは使わないので他は削除 ( 置いてあっても問題ない ) 名称を 0 0.org に変更 (0 のままだと snappyhexmesh でエラー ) 2013/8/10 22

23 BlockMesh 作成基本領域の 6 面体の設定とメッシュ分割 2013/8/10 23

24 blockmeshdict の格納場所 Cylinder フォルダの中身 blockmeshdict constant/polymesh の中に基本の 6 面体メッシュ作成を制御する blockmeshdict がある 2013/8/10 24

25 基本領域設定 (blockmeshdict) 節点の座標 (x y z) 6 面体指定変更箇所は赤枠単位変換 1m 0.01m // の後はコメント 2 次元問題とするため z 方向厚を 1mm と小さくする x,y,z 方向の分割数 0 z y 3 x 層流ではRe 9/4 個のセル数が目安単位長さあたりR 3/4! 20 3/4 =9.5なので円柱直径 10mm を10 分割する 1 セル 1mm 角 =1800 cells /8/10 25

26 基本境界 (Patch) 設定 (blockmeshdict) 7 upanddown 6 upstream 4 5 frontandback 3 0 frontandback 1 z y upanddown 2 downstream x 2 次元問題の時計算しない面 2013/8/10 26

27 blockmesh 作成コマンド実行フォルダへの移動 $cd cylinder blockmesh コマンド実行 $blockmesh 2013/8/10 27

28 blockmesh 実行後変更や追加されたファイル 2013/8/10 28

29 ParaFoam の起動 (BlockMesh の確認 ) ハイライトクリック ParaFoamの起動 Applyをクリック $parafoam 2013/8/10 29

30 メッシュの確認 Surface With Edges を指定 2013/8/10 30

31 Patch の確認 Patch 名が表示チェック 2013/8/10 31

32 snappyhexmesh 作成円柱周りのメッシュ作成 2013/8/10 32

33 ベースファイルの取得 tutorials/incompressible/pisofoam/les/motorbike/motorbike/system から以下の 3 ファイルをコピー 1snappyHexMeshDict (snappyhexmesh の制御ファイル ) 2decomposeDict ( 並列計算で使用これがないと snappy でエラー ) 3forceCoeffs(Cd 値計算で使用ついでにコピー ) 2013/8/10 33

34 snappyhexdict の変更 1 段階 ( ): 表面に適合したメッシュ除去 2 段階 ( ): 表面メッシュのスナップ ( 平滑化 ) 3 段階 ( ); 表面へのレイヤー挿入円柱ファイルの名称ファイルタイプ :stl 表面の名前指定 :cylinder 2013/8/10 34

35 snappyhexmeshdict の変更細分化したい表面の名称細分化の程度 ( 最小最大 ) 数が大きくなるほど細分化最小レベルは全ての表面に適用最大レベルは resolvefeatureangel 以上の角度の表面に適用細分化したい Box 領域の細分化今回使わないのでコメントアウト Snappy を適用する領域を指定今回は円柱の外部の座標を選ぶセルの表面の座標は好ましくない但し今回の点はセル表面であるが認識できた 2013/8/10 35

36 snappyhexmesh 実行 snappyhexmesh 実行 $snappyhexmesh 2013/8/10 36

37 snappyhexmesh 実行 snappy 成功 cell 数はディレクトリが生成 ( 後に説明するcotrolDictの時間ディレクトリ作成間隔が0.005secの場合 ) 2013/8/10 37

38 メッシュの確認 ( 第 1 段階 ) 次の時間に進める時間 $parafoamを起動円柱部分のセルが除去されている表面は階段状 2013/8/10 38

39 メッシュの確認 ( 第 2 段階 ) 時間 0.01 円柱表面の平滑化 2013/8/10 39

40 メッシュ確認円柱側面は1 層だけ切れている (2 次元問題 OK) 円柱側面のPatchはcylinder_zone0の名前が付けられている 2013/8/10 40

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

42 0 ディクショナリ作成初期値境界値の設定 2013/8/10 42

43 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

44 p=(p/ρ) 0.org/p 次元 =(m 2 /s 2 ) 初期値全域 0 p=0 p=0 p=0 2 次元問題条件 2013/8/10 あまり適切ではないが選択できる条件では 44 これが最もよい ( 付録参照 )

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

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

47 物性値と解析モデル turbulenceproperties ファイル ( 変更なし ) RASProperties ファイル乱流モデルの選択 RASmodel (Reynolds 平均モデル ) LESmodel (Large Eddy Simulation モデル ) tranportproperties ファイル ( 変更なし ) 動粘性係数 [m 2 /s] 輸送係数モデル選択 Newtonian ( ニュートン流体 ) CrossPowerLawCoefss ( 非ニュートン流体 ) など層流 2013/8/10 47

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

controldict 設定計算ソルバー starttime の設定時間から計算開始計算開始時間 endtime の設定時間に計算終了計算終了時間 ( 試行錯誤で決定 ) 計算時間ステップ ( クーラン数で決定 ) witeinterval ごとに時刻ディクショナリ生成 1000 ステップ 1

49 controldict 設定計算ソルバー starttime の設定時間から計算開始計算開始時間 endtime の設定時間に計算終了計算終了時間 ( 試行錯誤で決定 ) 計算時間ステップ ( クーラン数で決定 ) witeinterval ごとに時刻ディクショナリ生成 1000 ステップ 1 秒ごと時刻ディレクトリ最大数 : 制限なしデータファイルフォーマット選択 writeformatでの有効桁数データ圧縮指定 timeformatのフォーマット選択 timeformat の設定で用いる指数各ディクショナリの yes/no スイッチ物体に係る力の係数計算用の関数 2013/8/10 49

50 forcecoeffs 円柱 Patch 密度 [kg/m 3 ] 揚力方向 y 抗力方向 x ピッチモーメントの回転軸 z 上流での流速 [m/s] 物体中心代表寸法 [m] 射影面積 [m 2 ] 2013/8/10 50

51 fvscheme ( 変更なし ) 2013/8/10 51

52 fvscheme ( 変更なし ) 2013/8/10 52

53 Fvsolution( 変更なし ) 2013/8/10 53

54 STL ファイルの格納 2013/8/10 54

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

56 計算実行 2013/8/10 56

57 pisofoam の実行 $pyfoamplotrunner.py pisofoam pyfoamplotrunner.py で残差等がモニターできる 2013/8/10 57

58 計算終了時の画面計算時間抗力係数 ( 実験値は 2.0 程度 ) 揚力係数ピッチモーメント係数 2013/8/10 58

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

60 parafoam スケールをクリック U を選択時間経過を観察図の拡大マウス中軸を回す 2013/8/10 60

61 流線作図 StreamTracer クリッククリック 2013/8/10 61

62 流線作図 5 クリック 1Line Source を選択を入力 2Y Axis を選択 450 を入力 2013/8/10 62

63 流線作図 2013/8/10 63

64 速度分布と流線の合図目玉アイコンを ON Opacity( 透明度 ) 双子渦もそれらしい 2013/8/10 64

65 まとめ pisofoam を用いて 2 次元円柱周り流れの解析の演習を実施解析の流れは以下のとおり 1 円柱 STL ファイル作成 (FreeCAD) 2blockMesh ファイル作成実行 3snappyHexMesh ファイル作成実行 4 各種ファイル ( 境界物性制御 ) の作成 5pisoFoam の実行 6paraFoam を使った流線図の作成 2013/8/10 65

66 付録 Ⅰ 境界条件メッシュ影響 2013/8/10 66

67 境界条件 p の影響 Patch U p(case1) p(case2) p(case3) p(case4) Mesh upstream (2,0,0) 0 0 zerogradient zerogradient downstream (2,0,0) upanddown (2,0,0) 0 0 zerogradient zerogradient cylinder_zone0 (0,0,0) zerogradient fixedvalue; value $internalfield fixedvalue; value $internalfield zerogradient frontandback empty Cd 双子渦 - 計算時間 Cd 値は実験 (Cd=2.0) とかなり異なる特に円筒壁には $internalfieldを使うのは不適 2013/8/10 67

68 境界条件 p の影響 case1 case2 case3 case4 双子渦の様子が実験と異なる 2013/8/10 68

69 計算領域の影響 Patch U p(case1) p(case5) p(case6) Δ x-δ y - 6-3cm 24-12cm 28-16cm Base Mesh 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 双子渦 - 計算時間 ,404 2,054 4 倍程度領域が広くなるとCdが収束してくる但し計算時間は10 倍 2013/8/10 69

70 計算領域の影響 (case5) 2013/8/10 70

71 付録 Ⅱ Re=200 の場合 2013/8/10 71

72 領域とメッシュの影響項目 case1 case2 case3 Δ x-δ y 6-3cm 24-12cm 24-12cm Base Mesh Cells ,884 Cd St 計算時間 (sec) 129 (2.2min) 2,729 (45.5min) 19,987 (5.5hrs) 実験では Cd=1.5 St=0.19 なので Case3 では程々の精度 2013/8/10 72

73 Cd と Cl の履歴 (case3) 2013/8/10 73

74 流れの様子 (case3) 2013/8/10 74

75 付録 Ⅲ extrudemesh を使った 2D メッシュの作成 2013/8/10 75

76 概要オープン CAE 富山 / 第 10 回での富山大学の中川先生の資料を参考オープン CAE 勉強会 - 富山 /#_99 snappyhexmesh で 3 次元メッシュを作成 extrudemesh で Patch 表面メッシュだけを押し出して 2 次元メッシュとする演習での方法では円柱周りメッシュが 1 層にできない場合があるが本方法では必ず 1 層にできる本方法で領域のメッシュを部分的に詳細にすることが可能で計算時間の短縮になる 2013/8/10 76

77 blockmeshdict の修正分解 2013/8/10 77

78 snappyhexmeshdict の修正 1 領域名称 < 詳細メッシュ領域の指定 > 2 つのボックスを指定 refinementbox1 refinementbox2 ( 名前は任意 ) 指定図形タイプ : ボックスボックスを規定できる2 点を指定領域名称指定図形タイプ : ボックスボックスを規定できる 2 点を指定 2013/8/10 78

79 snappyhexmeshdict の修正 2 < 詳細メッシュ領域の分割法指定 > 領域名称ボックス内部 / 外部を指定領域名称細分化レベル levels( 距離レベル ) ボックス内部 / 外部を指定細分化レベル levels( 距離レベル ) refinementbox1 の細分化レベルが 1 に対し refinementbox2 はレベル 2 でより細かく細分化される 2013/8/10 79

80 extrudemeshdict の作成 pimpledymfoam/wingmotion/wingmotion2d_simplef oam/system/extrudemeshdict を /system にコピー現在のディレクトリでソースケースを参照押し出す Patch: front 押し出された Patch: back 押し出す層数 2013/8/10 80

81 U と p の修正 frontandback の分離 2013/8/10 81

82 snappyhexmesh 後のメッシュ $blockmesh $snappyhexmesh -overwrite 細かい領域だけ 5 層 2013/8/10 82

83 extrudemesh の実行後のメッシュ $extrudemesh 一層だけ 2013/8/10 83

84 Cd と St の比較項目今回等間隔メッシュ Δ x-δ y 24-12cm 24-12cm Cells 38, ,884 Cd St 計算時間 (sec) 1,880 (31.3min) 19,987 (5.5hrs) 2013/8/10 84

85 流線等間隔メッシュ今回メッシュ 2013/8/10 85

この講習会では, ユーザーガイドを参照しながら, 作業を進めますユーザーガイドは下記サイトで入手可能です OpenFOAM Documentation( オリジナル )

この講習会では, ユーザーガイドを参照しながら, 作業を進めますユーザーガイドは下記サイトで入手可能です OpenFOAM Documentation( オリジナル ) はじめての OpenFOAM その 2 富県学中川慎二オープンCAE 勉強会 @ 富 2014 年 1 月 25 日 Disclaimer: OPENFOAM is a registered trade mark of OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM and OpenCFD

PisoFoamによる 2次元円柱周りの流れの解析

PisoFoamによる　　　　　　　　　2次元円柱周りの流れの解析