こんにちは（@t_kun_kamakiri）
本記事ではOpenFOAMのメッシュ生成ユーティリティであるblockMesh+snappyHexMeshを使って円柱内部のメッシュ生成を行いたいと思います。
最終的には以下のようになります。
snappyHexMeshを使って円筒内部のメッシュ生成を行う
以前の記事でも解説したようにsnappyHexMeshでの境界層はなかなかきれいに生成することができません。
次回に解説するcfMeshの方が境界層はきれいに生成することができます。
ただ、cfMeshよりsnappyHexMeshの方が設定項目が多く融通が利きますが慣れが必要です。未だに使うかは悩ましいところですね。
cfMeshの境界層についてはこちらをご参考ください。
円筒内部の流れは理論的にもよく理解されている基本的な流れです。
- レイノルズ数の大きさによる層流・乱流遷移
- 層流と乱流の助走区間の違い
- 層流と乱流の速度分布の違い
これらの内容は流体力学の参考書にもある話です。
流体のシミュレーションをしていると何となく結果はでるけども「あってるのどうなの？」と思うこともしばしばあります。
実験とシミュレーションの結果を比較するというのもひとつの手ではありますし、流体力学の中でも条件を限定すると解析解を導くことができ、シミュレーションとの比較を行うことができます。
近似的に解析解が出せるなら、シミュレーションの結果と比較して妥当性を検証することは重要なことです。
今回はシミュレーションが正しく行えているのかを確認するために、円筒のモデルを作成して流体解析を行い、解析解とシミュレーションとの比較を行うことにします。
今回はメッシュだけを作ります。
【使用環境】
OpenFOAM-v2012(WSL2)
メッシュ生成の流れ
snappyHexMeshでのメッシュ生成には以下の手順が必要です。
- 形状ファイルを作成(FreeCAD)
- ベースメッシュの生成(blockMesh)
- 特徴線の抽出(surfaceFeatureExtract)
- 形状まわりのメッシュ生成(snappyHexMesh)
かっこ内は手法を書きましたが、オープンソースを使ったオーソドックスな方法を取っているに過ぎません。
形状ファイルはFusion360を使っても良いし、salomeを使っても良いです。
メッシュ生成までもsalomeで行っても良いです。
FreeCADで円柱モデルを作成
FreeCAD(0.20)でΦ10mm×1200mmの円柱モデルを作成します。
snappyHexMeshでの形状の認識のために形状stlファイルが必要であるためです。
また、FreeCADにはソリッドを面分割して面に名前を付ける機能があります。
ここで面に名前を付けておくことと便利です。面に名前を付けることで、境界条件の設定やsnappyHexMeshでの面まわりのメッシュの再分割などにも使えます。
面に名前をつけたstlファイルを結合してひとつのファイルにまとめるFreeCADのマクロは以下の記事をご参考ください。
出力したstlファイルは「constant/triSurface」に保存します。
FreeCADでの「mm単位」の表記は実は数値しか拾っていないので、実際流体解析で使用する単位系がm単位であればスケール変換する必要があります。
スケール変換は以下のコマンドで簡単に行うことができます。
1 | surfaceConvert -scale 0.001 constant/triSurface/model.stl constant/triSurface/model_m.stl |
今回使用するのはm単位系になった「model_m.stl」です。
チュートリアルをコピー
適当なチュートリアルからコピーしてきます。
1 | cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily pipe_snappyHexMesh |
ここで、環境変数「$FOAM_TUTORIALS」が「/opt/OpenFOAM/OpenFOAM-v2012/tutorials」になっています。
フォルダを移動します。
1 | cd pipe_snappyHexMesh |
ベースメッシュの生成(blockMesh)
まずはベースメッシュから生成します。
ベースメッシュはOpenFOAMのユーティリティであるblockMeshを使います。blockMeshの設定ファイルは「system/blockMeshDict」で行います。
PraViewでstlファイルの寸法を確認しながら設定しましょう。
system/blockMeshDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object blockMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // scale 0.001; // Divisions in x/y/z directions. Can be unequal. nx 10; ny 10; nz 500; Xmin -11.1; Xmax 11.1; Ymin -11.1; Ymax 11.1; Zmin -1.1; Zmax 1201.1; geometry { } vertices ( ($Xmin $Ymin $Zmin) //0 ($Xmax $Ymin $Zmin) //1 ($Xmax $Ymax $Zmin) //2 ($Xmin $Ymax $Zmin) //3 ($Xmin $Ymin $Zmax) //4 ($Xmax $Ymin $Zmax) //5 ($Xmax $Ymax $Zmax) //6 ($Xmin $Ymax $Zmax) //7 ); blocks ( hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1) ); edges ( ); faces ( ); boundary ( ); // ************************************************************************* // |
面の名前はFreeCADで付けましたので、boundaryで面の名前を指定する必要はありません。また、境界のタイプはsnappyHexMeshで行いますので境界のタイプの指定も不要です。
あくまでblockMeshはsnappyHexMeshのベースのメッシュのためだけに行うものです。
メッシュのある程度のサイズはblockMeshで決まってしまうため、メッシュサイズは少し注意を払う必要があります。
特徴線の抽出(surfaceFeatureExtract)
今回のような単純な形状の場合は特徴線の抽出の恩恵を受けることができないですが、形状が複雑な場合には特徴線を抽出しておくとsnappyHexMeshの設定で特徴線まわりでメッシュの再分割をすることができるのでメッシュ生成が楽になります。
特徴線の設定は「system/surfaceFeatureExtractDict」で行います。
しかし、コピーしてきたチュートリアルにはsurfaceFeatureExtractDictがないため、適当なチュートリアルから持ってくる必要があります。
どのチュートリアルにsurfaceFeatureExtractDictがあるのかを知っていれば良いですが、知らない場合が多いので以下のコマンドで探すことにします。
1 | find $FOAM_TUTORIALS -name "surfaceFeatureExtractDict" |
するといくつか候補が出てきます。
1 2 3 4 5 6 | /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/surfaceFeatureExtractDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/motorBike/system/surfaceFeatureExtractDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/mixerVesselAMI/system/surfaceFeatureExtractDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/DTCHullMoving/system/surfaceFeatureExtractDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/DTCHull/system/surfaceFeatureExtractDict (以下省略) |
適当にコピーしてsystemに保存します。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/simpleFoam/motorBike/system/surfaceFeatureExtractDict system/ |
system/surfaceFeatureExtractDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object surfaceFeatureExtractDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // model_m.stl { // How to obtain raw features (extractFromFile || extractFromSurface) extractionMethod extractFromSurface; // Mark edges whose adjacent surface normals are at an angle less // than includedAngle as features // - 0 : selects no edges // - 180: selects all edges includedAngle 150; subsetFeatures { // Keep nonManifold edges (edges with >2 connected faces) nonManifoldEdges no; // Keep open edges (edges with 1 connected face) openEdges yes; } // Write options // Write features to obj format for postprocessing writeObj yes; } // ************************************************************************* // |
設定はそれほど難しくないですね。
どの形状まわりで特徴線を抽出するのか「model_m.stl」を指定して、特徴線を判定する角度「includedAngle 150;」を指定するだけです。
では、以下のコマンドで特徴線を生成します。
1 | surfaceFeatureExtract |
コマンド実行後に「constant/triSurface」の中に「model_m.eMesh」というファイルができます。
こちらが特徴線のファイルになりますが、ParaViewで直接読み込むことができないので以下のコマンドでParaViewで読み込める形式に変換します。
1 | surfaceFeatureConvert constant/triSurface/model_m.eMesh constant/triSurface/model_m.obj |
ParaViewで読み込むと以下のように円柱の円形の部分だけが表示されます。
形状まわりのメッシュ生成(snappyHexMesh)
形状ファイル、ベースメッシュ、特徴線ができたらsnappyHexMeshで形状まわりのメッシュを生成します。
snappyHexMeshは「system/snappyHexMeshDict」で設定を行いますが、今回コピーしてきたチュートリアルにはなかったので、また適当なチュートリアルから持ってくる必要があります。
以下のコマンドでsnappyHexMeshDictが使われているチュートリアルを探し、
1 | find $FOAM_TUTORIALS -name "snappyHexMeshDict" |
いくつか候補が出てきます。
1 2 3 4 5 6 7 | /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/snappyHexMeshDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/snappyHexMeshDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/motorBike/system/snappyHexMeshDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/mixerVesselAMI/system/snappyHexMeshDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/DTCHullMoving/system/snappyHexMeshDict /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/DTCHull/system/snappyHexMeshDict (以下省略) |
適当なチュートリアルからsystemにコピーします。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict system/ |
あともう一つ、snappyHexMeshDictにメッシュ品質（meshQualityDict）に関する記述がありますが、その設定ファイルも持ってくる必要があります。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/simpleFoam/motorBike/system/meshQualityDict system/ |
では、snappyHexMeshDictの設定に移ります。
system/snappyHexMeshDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object snappyHexMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Which of the steps to run castellatedMesh true; snap true; addLayers true; geometry { model_m.stl { type triSurfaceMesh; name model_m; regions { inlet { name inlet; } outlet { name outlet; } sideWall { name sideWall; } } } // refinementBox // { // type box; // min (-1.0 -0.7 0.0); // max ( 8.0 0.7 2.5); // } } // Settings for the castellatedMesh generation. castellatedMeshControls { features ( { file "model_m.eMesh"; level 0; } ); refinementSurfaces { model_m { level (0 0); regions { inlet { level (0 0); patchInfo { type patch; // inGroups (motorBikeGroup); } } outlet { level (0 0); patchInfo { type wall; } } sideWall { level (1 2); patchInfo { type wall; } } } } } refinementRegions { // refinementBox // { // mode inside; // levels ((1E15 4)); // } } maxLocalCells 100000000; maxGlobalCells 100000000; minRefinementCells 10; maxLoadUnbalance 0.10; nCellsBetweenLevels 3; resolveFeatureAngle 30; locationInMesh (0 0 0.1); allowFreeStandingZoneFaces true; } // Settings for the snapping. snapControls { nSmoothPatch 3; tolerance 2.0; nSolveIter 30; nRelaxIter 5; nFeatureSnapIter 10; implicitFeatureSnap false; explicitFeatureSnap true; multiRegionFeatureSnap false; } // Settings for the layer addition. addLayersControls { layers { //"(lowerWall|motorBike).*" inlet { nSurfaceLayers 3; } outlet { nSurfaceLayers 3; } sideWall { nSurfaceLayers 3; } } relativeSizes true; expansionRatio 1.0; finalLayerThickness 0.3; minThickness 0.03; nGrow 0; // featureAngle 60; featureAngle 130; slipFeatureAngle 30; nRelaxIter 3; nSmoothSurfaceNormals 1; nSmoothNormals 3; nSmoothThickness 10; maxFaceThicknessRatio 0.5; maxThicknessToMedialRatio 0.3; minMedialAxisAngle 90; nBufferCellsNoExtrude 0; nLayerIter 50; } // Generic mesh quality settings. At any undoable phase these determine // where to undo. meshQualityControls { #include "meshQualityDict" //- Number of error distribution iterations nSmoothScale 4; //- Amount to scale back displacement at error points errorReduction 0.75; } // Write flags writeFlags ( scalarLevels layerSets layerFields // write volScalarField for layer coverage ); mergeTolerance 1e-6; // ************************************************************************* // |
設定が多くてややこしいですが以下の部分だけを押さえておけば良いでしょう。
- geometry：形状ファイルの認識（面の名前も指定）
- castellatedMeshControls.features：特徴線まわりの再分割（再分割しないのでレベル0)
- castellatedMeshControls.refinementSurfaces：指定した面のまわりで再分割
- castellatedMeshControls.refinementRegions：指定した領域で再分割(今回は設定なし）
- addLayersControls.layers：指定した面で境界層を設定
では、snappyHexMeshを実行してメッシュを作ります。
1 | snappyHexMesh -overwrite |
オプションで「-overwrite」を付けるとメッシュ生成の途中経過を上書きしてくれます。
snappyHexMeshは3段階でメッシュを生成するため、オプションなしで実行するとメッシュ生成の途中までフォルダ分けして作られます（「1」「2」「3」のように・・・「3」が最終状態なので「3/polyMesh」を「constant」に移動すれば良いですが、面倒なので途中経過が必要ない方は「-overwrite」を付けて実行した方が良いでしょう。
メッシュ生成には数分時間がかかります。
ParaViewで確認するとこのようになっています。
断面を切るとやはり境界層がきれいではないですね。snappyHexMeshの仕様の問題で、仕方ないですね・・・・
メッシュ品質の確認
メッシュ品質も確認しておきましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 | Mesh stats points: 347747 faces: 981408 internal faces: 933016 cells: 317088 faces per cell: 6.03752 boundary patches: 3 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 283008 prisms: 24096 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 9984 Breakdown of polyhedra by number of faces: faces number of cells 6 3976 8 32 12 5952 15 24 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology sideWall 47904 47952 ok (non-closed singly connected) outlet 244 261 ok (non-closed singly connected) inlet 244 261 ok (non-closed singly connected) Checking faceZone topology for multiply connected surfaces... No faceZones found. Checking basic cellZone addressing... No cellZones found. Checking geometry... Overall domain bounding box (-0.00999935 -0.00999906 0) (0.00999935 0.01 1.2) Mesh has 3 geometric (non-empty/wedge) directions (1 1 1) Mesh has 3 solution (non-empty) directions (1 1 1) Boundary openness (-5.38986e-15 3.37246e-14 2.21912e-18) OK. Max cell openness = 3.6806e-16 OK. Max aspect ratio = 7.86712 OK. Minimum face area = 8.22387e-08. Maximum face area = 5.58801e-06. Face area magnitudes OK. Min volume = 4.66207e-11. Max volume = 1.24305e-08. Total volume = 0.000375705. Cell volumes OK. Mesh non-orthogonality Max: 43.0107 average: 7.37739 Non-orthogonality check OK. Face pyramids OK. Max skewness = 3.43072 OK. Coupled point location match (average 0) OK. Mesh OK. |
「faces: 981408」なので面の数が90万くらいあるということですね。
個人のPCで定常の流体解析をすると並列計算（手元のPCは6並列）でも数時間かかる感じですね。
以下の部分は最低確認するようにしましょう。
- Max aspect ratio：アスペクト比（縦横の比率）→境界層で大きくなりやすいので大きすぎる場合は確認
- non-orthogonality：直交性→70°以上であれば修正
- Max skewness：歪度→4以上だと修正
メッシュ品質の基準は以下のファイルにも記載がありますで確認してみましょう。
1 | $WM_PROJECT_DIR/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/ |
コマンドでファイルを探す場合は、
1 | find $FOAM_SRC -name "primitiveMeshCheck*" |
以下のように出力されます。
1 2 3 4 5 6 7 | /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude/primitiveMeshCheckEdgeLength.C /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude/primitiveMeshCheckPointNearness.C /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude/primitiveMeshCheck.C /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckEdgeLength.C /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckPointNearness.C /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C |
primitiveMeshCheck.Cの中身を確認しましょう。
1 | vi $FOAM_SRC/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C |
該当箇所は「/primitiveMesh」と打って「n(下へ）」「shift + n(上へ）」で探すことができます。
viを終了するには
ESCボタンを押して「:q」か「:q!」で終了することができます。
1 2 3 4 5 | Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6; Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000; Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4; Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6; |
シミュレーションを行う場合、低品質なセルや面によって生じる数値誤差を低減できる数値スキームを選択します。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
Follow @t_kun_kamakiri