OpenFAOM_pisoFoamによる2次元円柱周りの解析

Transcript

1. OpenFAOM_pisoFoamによる2次元
   円柱周りの解析
   2017年8月14日
   【参考資料】
   http://opencae.gifu-
   nct.ac.jp/pukiwiki/index.php?plugin=attach&refer=%C2%E8%A3%B2%A3%B4%B2%F3%CA%D9%B6%AF
   %B2%F1%A1%A7H250810&openfile=2Dflow_around_cylinder_with_pisoFoam_rev.pdf
   http://www.mech.iwate-u.ac.jp/~hirose/ockitatohoku/ref/kasaisensei-text-1.pdf
   

   View Slide

2. 目次
   1. 理論
   2. 検証項目
   3. OpenFOAMによる作成手順
   4. 結果
   5. 考察
   

   View Slide

3. 3
    2次元円柱回りの流れ
   基礎方程式（運動量保存）
   𝜕𝒖
   𝜕𝑡
   + 𝒖 ∙ 𝜵 𝒖 = −
   1
   𝜌
   𝜵𝑃 + 𝜈𝜵2𝒖
   𝜵 ∙ 𝒖 = 𝟎
   非圧縮の条件式(質量流保存より)
   代表的な長さ𝐷、速度𝑈、時
   間𝜏を用いて無次元化を行う 𝜕𝒖∗
   𝜕𝜏
   + 𝒖∗ ∙ 𝜵∗ 𝒖∗ = −𝜵∗𝑃∗ +
   1
   𝑅𝑒
   𝜵∗2𝒖∗
   レイノルズ数：𝑅𝑒 =
   𝑈𝐷
   𝜈
   無次元化することにより、基礎方程式の
   パラメータはレイノルズ数のみとなる
   レイノルズ数が同じであれば、流れ場は
   同じ振る舞いをする。(力学的相似則)
   𝜵∗ ∙ 𝒖∗ = 𝟎
   非圧縮の条件式(質量流保存より)
   基礎方程式（運動量保存）
   レイノルズ数Reによる流れ場の変化
   ⚫ 非圧縮性流体の支配方程式
   ⚫ レイノルズ数と流れ場
   

   View Slide

4.  OpenFOAMによる作成手順
   円柱stlファイル作成（FreeCAD)
   メッシュ作成Ⅰ（blockMesh)
   メッシュ作成Ⅱ（snappyHexMesh)
   境界条件、計算条件の設定
   計算実行
   結果処理
   プリ処理
   ソルバー
   ポスト処理
   

   View Slide

5. 円柱stlファイル作成（FreeCAD)
   ①「FreeCAD」で円柱を作成
   半径：5.0mm
   高さ：10.0mm
   ②メッシュの作成
   「メッシュ（M)」→「シェイプ
   からメッシュを作成する」
   Meshing Tolerance
   ここでは「標準」を選択し、デ
   フォルト設定のままでOK
   

   View Slide

6. ③「メッシュ(M)」→「メッシュ
   をエクスポート」を選択
   ファイル名：cylinder.stl
   ファイル形式：Binary STL
   (1)
   (2)
   ④単位変換
   FreeCADで10mmと設定したが、
   mmは認識されていない。
   ※10ｍと認識されている。
   10/1000=0.01m
   と単位変換するため、単位変換コ
   マンドを実行
   $surfaceConvert -scale 0.001 cylinder.stl cylinder-s.stl
   もしくは、
   $surfaceTransformPoints cylinder.stl cylinder-s.stl -scale ’(0.001 0.001 0.001)’
   

   View Slide

7. ①ベースファイルの取得
   ※OpenFOAMでは現象にあったソルバーをチュートリアルからコピー
   して計算条件の設定を変更して解析を実行する。
   メッシュ作成Ⅰ（blockMesh)
   ベースとなるチュートリアル保存先
   ホーム/Desktop/OpenFOAM-4x/tutorials/incompressible/pisoFoam/ras/cavity
   これを作業フォルダにコピー
   ファイル構造
   0
   constant
   system
   時刻ディクショナリ
   U,Pなどの物理量の初期値、境界条件
   constantディクショナリ
   transportProperties：輸送物性値
   turbulenceProperties：乱流モデル（RANS、LES）指定
   systemディクショナリ
   blockMeshDict：ベースメッシュ
   controlDict：ソルバー、時間設定（終了時間、タイムステップ）
   fvSchemes：有限体積法の離散化方法の指定
   fvSolution：有限体積法の方程式解放と収束判定の指定
   

   View Slide

8. ②「0」ディレクトリの名前を変更
   ここでは「0.org」に変更した。
   ※名前が「0」のままだと
   snappyHexMeshでエラーが出る。
   「0」ディレクトリの中身
   U:速度
   P：圧力
   epsilon：乱流エネルギー散逸率
   k：乱流エネルギー
   nuTilda：渦粘性：Spalart-Allmarasモデル
   （1方程式モデル）
   nut：渦粘性νt
   omega：ω
   

   View Slide

9. ③基本領域の設定
   ファイル：System/blockMeshDict
   0 1
   2
   6
   7
   4
   3
   5
   乱流におけるセル分割数
   𝑳
   𝜼
   ~𝑹𝒆
   𝟑
   𝟒→Re=20とすると20
   3
   4 = 9.5
   円柱10mmを10分割
   ※分割数はエネルギー供給率と散逸率が等し
   いとして、次元解析より見積もる
   x
   y
   z
   単位変換 1m→0.01m
   30mm
   1mm
   60mm
   各ボックスの座標を指定
   //はコメント
   メッシュの勾配
   （１１１）：均一
   (0.5,1,1)：x方向にx方
   向の分割数が小さくなる
   x,y,z方向の分割数
   ６面体指定
   

   View Slide

10. ④境界条件の設定
    辺の形状を定義
    ここでは使用しない
    境界定義
    0 1
    2
    6
    7
    4
    3
    5
    upANDdown
    upANDdown
    downstream
    upstream
    frontANDback
    face番号
    empty
    2次元
    問題の
    時には
    計算し
    ない面
    face名
    境界のタイプ
    patch パッチ
    wall 壁
    symmetryPlane 対称面
    cyclic 周期境界
    cyclicAMI 不整合周期境界
    wedge 2次元軸対称
    empty 2次元
    結合する境界面の対を指定
    例）境界名patch1,patch2の時、(patch1, patch2)とする
    

    View Slide

11. ⑤blockMeshの作成
    $cd 作業フォルダ名
    $blocMesh
    作業フォルダでblockMesh
    コマンドを実行するとこの階層にあ
    る「blockMeshDict」ファイルを実
    行してくれる。
    blockMeshを実行すると、
    constant/polyMesh
    にnode,face,boundaryなどの情報ファイル
    が作成される。
    

    View Slide

12. ⑥blockMeshの確認
    ＄paraFoam （paraFoamの起動）
    Patch Names
    Patch名表示
    「Surface With Edges」
    メッシュ表示
    

    View Slide

13. メッシュ作成Ⅱ（snappyHexMesh)
    ①ベースファイルから以下の3つの
    ファイルを作業フォルダにコピー
    ベースファイル
    tutorials/incompressible/pis
    oFoam/les/motorBike/moto
    rBike/system
    (1)snappyHexMeshDict
    snappyHexMeshの制御ファイル
    (2)decomposeParDict
    並列計算
    (3)forceCoeffs
    抗力係数Cdを計算
    

    View Slide

14. snappyHexMeshの仕組み
    次のような手順でメッシュを作成する。
    STLなどのサーフェスを用意。 blockMeshでベースとなる
    メッシュを作成。
    snappyHexMeshで境界近傍
    を六面体で細分化する。図で
    は一部分だけ実施した。
    サーフェイス全体に対
    して行う。
    領域外になる部分のセ
    ルを取り除く。
    領域外に飛び出た部分を
    サーフェイスに合わせる
    ("snap")
    

    View Slide

15. 1段階（○）：表面に適合したメッシュ除去
    2段階（○）：表面メッシュのスナップ（平滑化）
    3段階（×）：表面へのレイヤー挿入
    円柱のファイル名
    ②snappyHexMeshDictの編集
    ファイルタイプ：stl
    表面の名前指定：cylinder
    細分化したい表面の名称 細分化の程度（最小、最大）
    0⇒1⇒2⇒・・・大きくなるほど細分化
    • 最小レベルは全ての表面に適用
    • 最大レベルはresolveFeatureAngle以
    上の角度の表面に適応
    System/ snappyHexMeshDict
    

    View Slide

16. 細分化したいBox領域の細分化
    今回は使わないためコメントアウト
    Snappyを適用する領域を指定。今回は円柱
    の外部の座標を選ぶ。セルの表面の座標は好
    ましくない。但し、今回の点はセル表面であ
    るが、認識できた。
    0 1
    2
    6
    7
    4
    3
    5
    0(-0.002,-1.5,0)
    円柱表面の座標
    

    View Slide

17. ③「constant」フォルダに新規でtriSurfaceという
    名前のディレクトリを作成し、円柱形状の
    cylinder-s.stlファイルを置く
    ※この操作をしない場合は、$snappyHexMeshを
    実行すると下記のようなエラーが出される。
    ④snappyHexMeshの実行
    $snappyHexMesh
    メッシュ作成の計算がエラーなく実
    行された場合は左のよう書き出され
    終了する。
    

    View Slide

18. ⑤メッシュの確認
    $paraFoamでメッシュの確認
    第一段階
    円柱部分のセルが除去されている
    第二段階
    円柱表面の平滑化がされている
    時間 0.01
    時間 0.0
    円柱側面は1層だけ切れ
    ている。
    円柱側面のPatch名
    は”cylinder”
    

    View Slide

19. ⑥生成したディレクトリを削除
    最終時間の0.01を解析の初期状態とする。
    $rm –r 0.005
    $mv 0.01 0
    $mv 0.org/U 0
    $mv 0.org/P 0
    0.05を削除
    0.01のディレクトリの名前を0に変更
    0.orgディレクトリの中のU,Pを0ディレク
    トリに移す。（後で初期条件を設定）
    

    View Slide

20. 境界条件、計算条件の設定
    0 1
    2
    6
    7
    4
    3
    5
    upANDdown
    upANDdown
    downstream
    upstream
    frontANDback
    次元
    [kg m s K mol A cd]
    領域初期値 全域0.02m/s
    固定値 上流0.02m/s
    流出条件
    流入条件
    流出時
    流出条件
    流入条件
    流出時
    2次元問題
    固定値 粘着条件
    0/polyMesh/boundary 0/U
    ①速度条件の設定
    

    View Slide

21. 0/P
    ②圧力条件の設定
    次元[m2/s2]
    初期値 全域0
    p=0
    p=0
    p=0
    2次元問題条件
    𝜕𝑃
    𝜕𝑛
    = 0
    

    View Slide

22. ③物性値の選択
    動粘性係数
    次元[m2/s] 1.e-5
    輸送係数モデルの選択
    Newtonian ニュートン流体
    powerLaw
    非ニュートン流体
    BirdCarreau
    CrossPowerLaw
    HerschelBulkley
    ニュートン流体
    せん断力が速度勾配に比例
    𝜏 = 𝜇
    𝜕𝑢
    𝜕𝑦
    0/constant/transportProperties
    

    View Slide

23. ④流体モデルの選定
    simulationType モデルの種類
    laminar(層流)
    LES(LES乱流モデル) DeardorffDiffStress
    Smagorinsky
    SpalartAllmarasDDES
    SpalartAllmarasDES
    SpalartAllmarasIDDES
    WALE
    dynamicKEqn
    dynamicLagrangian
    kEqn
    kOmegaSSTDES
    RAS(RANS乱流モデル) LRR
    LamBremhorstKE
    LaunderSharmaKE
    LienCubicKE
    LienLeschziner
    RNGkEpsilon
    SSG
    ShihQuadraticKE
    SpalartAllmaras
    kEpsilon
    kOmega
    kOmegaSST
    kOmegaSSTSAS
    kkLOmega
    qZeta
    realizableKE
    v2f
    0/constant/turbulenceProperties
    Laminar(層流)の場
    合は設定の必要なし
    simulationTypeが
    「RAS」,「LES」の場合は右表を参照
    ※円柱回りの流れはレイノルズ数が約2300以上で乱流に遷移する。
    乱流現象の場合は乱流モデル（RANS、LES）もしくはDNSでの計算を
    行う必要がある。
    

    View Slide

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

    View Slide

25. 物体 寸法 抗力係数
    ⑥抗力の設定
    0/constant/forceCoeffs
    円柱のPatch
    揚力y
    抗力x
    密度[kg/m3]
    物体中心
    ピッチモーメント回転軸z
    上流の流速[m/s]
    物体寸法[m]
    射影面積[m2]
    𝑣
    密度𝜌
    𝑑
    𝑙
    𝑙
    𝑑
    =
    1
    2
    5
    10
    40
    ∞
    0.63
    0.68
    0.74
    0.82
    0.98
    1.20
    𝐶𝐷
    =
    𝐷
    1
    2 𝜌𝑣2𝐴
    𝐴 = 𝑑𝑙
    D:物体寸法
    A：射影面積
    

    View Slide

26. ⑦スキームの設定
    0/constant/fvSchemes
    

    View Slide

27. ⑧解法の設定
    0/constant/fvSolution
    

    View Slide

28. 計算実行
    その他の主なソルバーは以下の通り。
    ➢ basic : 基本
    laplacianFoam : 拡散方程式ソルバー
    potentialFoam : ポテンシャル流れソルバー
    scalarTransportFoam : スカラー輸送方程式ソルバー
    ➢ incompressible : 非圧縮性流体
    icoFoam : 非定常層流解析ソルバー
    simpleFoam : 定常乱流解析ソルバー (SIMPLE 法)
    pisoFoam : 非定常乱流解析ソルバー (PISO 法)
    pimpleFoam : 非定常乱流解析ソルバー (PIMPLE = PISO + SIMPLE 法)
    heatTransfer : 熱流動
    buoyantBoussinesqSimpleFoam : 定常熱流動解析ソルバー (Boussinesq 近似)
    buoyantBoussinesqPimpleFoam : 非定常熱流動解析ソルバー (Boussinesq 近似)
    buoyantSimpleFoam : 定常熱流動解析ソルバー
    buoyantPimpleFoam : 非定常熱流動解析ソルバー
    ➢ multiphase : 混相流
    interFoam : VOF 法による 2 相流解析ソルバー
    multiphaseInterFoam : VOF 法による多相流解析ソルバー
    ①解析実行コマンド
    $pisoFoam
    

    View Slide

29. 計算実行
    ①解析実行コマンド（残差の確認）
    $pyFoamPlotRunner.py pisoFoam
    logファイル
    抗力係数3.04426
    理論値≒2.1
    

    View Slide

30. 結果処理
    U(速度)の表示
    透明度
    「StreamTracer」
    をクリック
    流線表示
    流速コンター図
    High Resolution
    Line Source
    解像度の部分を設定
    ①「流速コンター図」「流線表示」
    paraview
    

    View Slide

31. 補足
    2.0m/s
    レイノルス数
    Re＝UL/ν
    =2.0*0.01/(1.0-5)
    =2000
    乱流遷移レイノルズ数
    2300
    速度コンター図
    流線
    速度ベクトル
    

    View Slide

32. 32
     2次元円柱回りの流れ
    基礎方程式（運動量保存）
    𝜕𝒖
    𝜕𝑡
    + 𝒖 ∙ 𝜵 𝒖 = −
    1
    𝜌
    𝜵𝑃 + 𝜈𝜵2𝒖
    𝜵 ∙ 𝒖 = 𝟎
    非圧縮の条件式(質量流保存より)
    代表的な長さ𝐷、速度𝑈、時
    間𝜏を用いて無次元化を行う 𝜕𝒖∗
    𝜕𝜏
    + 𝒖∗ ∙ 𝜵∗ 𝒖∗ = −𝜵∗𝑃∗ +
    1
    𝑅𝑒
    𝜵∗2𝒖∗
    レイノルズ数：𝑅𝑒 =
    𝑈𝐷
    𝜈
    無次元化することにより、基礎方程式の
    パラメータはレイノルズ数のみとなる
    レイノルズ数が同じであれば、流れ場は
    同じ振る舞いをする。(力学的相似則)
    𝜵∗ ∙ 𝒖∗ = 𝟎
    非圧縮の条件式(質量流保存より)
    基礎方程式（運動量保存）
    レイノルズ数Reによる流れ場の変化
    ⚫ 非圧縮性流体の支配方程式
    ⚫ レイノルズ数と流れ場
    

    View Slide

33. 円柱の直径L
    240mm
    計算領域の拡大
    10mm 30mm
    60mm
    4倍
    120mm
    細かくメッシュを切る
    •レイノルズ数：Re=(0.01*2.0)/(1.0-5)=2000
    •セルの分割数
    L/η～Re(3/4)⇒Re＝2000とすると、
    η=0.24/300=0.8mmが最小渦サイズ
    流速
    [m/s]
    60mm
    レイヤー層：5層
    

    View Slide

34. View Slide

35. 0 1
    2
    6
    7
    4
    3
    5
    upANDdown
    upANDdown
    downstream
    upstream
    frontANDback
    セルサイズ
    120/60=2mm
    blockMesh
    

    View Slide

36. セルサイズ
    120/60=2mm
    ↓
    2mm/4=0.5mm
    0.8mmの最小渦サイズができるの
    で、それより小さいセル要素を準備
    2分割ｗ2回繰り返
    したので、4層に
    なっている
    外側は1層
    z方向の層が1層ではないので、2次
    元解析の計算ができない。
    snappyHexMesh
    

    View Slide

37. snappyHexMesh
    境界層の挿入
    「true」に
    する
    addLayersControls
    {
    relativeSizes true;
    // Per final patch (so not geometry!) the layer information
    layers
    {
    “cylinder"
    {
    nSurfaceLayers 3;
    }
    }
    expansionRatio 1.1;
    finalLayerThickness 0.3;
    minThickness 0.1;
    nGrow 0;
    featureAngle 60;
    nRelaxIter 3;
    nSmoothSurfaceNormals 1;
    nSmoothNormals 3;
    nSmoothThickness 2;
    maxFaceThicknessRatio 0.5;
    maxThicknessToMedialRatio 0.3;
    minMedianAxisAngle 90;
    nBufferCellsNoExtrude 0;
    nLayerIter 50;
    }
    層の拡大率
    3層目のメッシュ
    サイズ
    3層の境界層が作成された。
    境界層の面の名前を指定
    

    View Slide

38. pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam
    /system/extrudeMeshDictを/system
    にコピー
    // What to extrude:
    // patch : from patch of another case ('sourceCase')
    // mesh : as above but with original case included
    // surface : from externally read surface
    constructFrom patch;
    sourceCase ".";
    sourcePatches (front);
    // If construct from patch: patch to use for back (can be same as sourcePatch)
    exposedPatchName back;
    // Flip surface normals before usage. Valid only for extrude from surface or
    // patch.
    flipNormals false;
    //- Linear extrusion in point-normal direction
    extrudeModel linearNormal;
    nLayers 1;
    expansionRatio 1.0;
    linearNormalCoeffs
    {
    thickness 0.001;
    }
    押し出す面
    この面まで押し出す
    1層だけにする
    $extrudeMesh
    1層の
    メッシュ
    になった
    

    View Slide

39. constant/transportProperties
    動粘性係数ν=1.0e-5
    constant/turbulenceProperties
    層流条件
    

    View Slide

40. dimensions [0 1 -1 0 0 0 0];
    internalField uniform (0.2 0 0);
    boundaryField
    {
    in
    {
    type fixedValue;
    value uniform (0.2 0 0);
    }
    out
    {
    type inletOutlet;
    inletValue uniform (0.2 0 0)
    value $internalField;
    }
    side
    {
    type inletOutlet;
    inletValue uniform (0.2 0 0)
    value $internalField;
    }
    front
    {
    type empty;
    }
    back
    {
    type empty;
    }
    cylinder
    {
    type fixedValue;
    value uniform (0 0 0);
    }
    }
    0/U
    壁：滑り無し
    流速規定
    フリースリップ
    or流速を指定 zeroGra
    dient
    or
    outletIn
    letで逆
    流防止
    inletOutletは、ドメインに流入するときディリクレ型、流出するときzeroGradient
    outletInletは、ドメインから流出するときディリクレ型、流入するときzeroGradient
    

    View Slide

41. dimensions [0 2 -2 0 0 0 0];
    internalField uniform 0;
    boundaryField
    {
    in
    {
    type fixedValue;
    value uniform 0;
    }
    out
    {
    type zeroGradient;
    }
    side
    {
    type fixedValue;
    value uniform 0;
    }
    front
    {
    type empty;
    }
    back
    {
    type empty;
    }
    cylinder
    {
    type zeroGradient;
    }
    }
    0/p
    圧力は指
    定しない
    圧力は指
    定しない
    圧力規定
    圧力規定
    

    View Slide

42. application pisoFoam;
    startFrom startTime;
    startTime 0;
    stopAt endTime;
    endTime 5;
    deltaT 0.0001;
    writeControl timeStep;
    writeInterval 1000;
    purgeWrite 0;
    writeFormat ascii;
    writePrecision 6;
    writeCompression off;
    timeFormat general;
    timePrecision 6;
    runTimeModifiable true;
    adjustTimeStep on;
    maxCo 0.9;
    functions
    {
    #include "forceCoeffs"
    }
    system/controlDict
    クーラン数指定
    

    View Slide

43. $tail -n 10 log
    $pyFoamPlotRunner.py pisoFoam >log &
    流速＝
    2.0m/s
    圧力だけ変動が大きい。
    おそらく、非圧縮条件の式を満たすよう
    に圧力を決めているので、圧力だけ残差
    の変動が大きくなる。
    

    View Slide

44. 結果処理
    抗力係数Cdが
    2.1に近づく
    

    View Slide

45. 途中から計算
    

    View Slide

46. APPENDIX
    

    View Slide

47. Valid searchableSurface types :
    11
    (
    closedTriSurfaceMesh
    distributedTriSurfaceMesh
    searchableBox
    searchableCylinder
    searchableDisk
    searchablePlane
    searchablePlate
    searchableSphere
    searchableSurfaceCollection
    searchableSurfaceWithGaps
    triSurfaceMesh
    )
    snappyHexMeshDict
    

    View Slide

48. $ transformPoints -scale "(1000 1000 1000)"
    

    View Slide