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

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

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

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

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

③「メッシュ(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)’

①ベースファイルの取得 ※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：有限体積法の方程式解放と収束判定の指定

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

③基本領域の設定 ファイル：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方向の分割数 ６面体指定

④境界条件の設定 辺の形状を定義 ここでは使用しない 境界定義 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)とする

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

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

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

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

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

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

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

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

⑥生成したディレクトリを削除 最終時間の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ディレク トリに移す。（後で初期条件を設定）

境界条件、計算条件の設定 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 ①速度条件の設定

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

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

④流体モデルの選定 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での計算を 行う必要がある。

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

物体 寸法 抗力係数 ⑥抗力の設定 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：射影面積

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

⑧解法の設定 0/constant/fvSolution

計算実行 その他の主なソルバーは以下の通り。 ➢ 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

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

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

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

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

円柱の直径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層

No content

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

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

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層の境界層が作成された。 境界層の面の名前を指定

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層の メッシュ になった

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

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

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 圧力は指 定しない 圧力は指 定しない 圧力規定 圧力規定

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 クーラン数指定

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

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

途中から計算

APPENDIX

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

$ transformPoints -scale "(1000 1000 1000)"