【OpenFOAM】snappyHexMeshで円柱周り流れのメッシュを作成する方法

#はじめに
この記事ではsnappyHexMeshユーティリティを使ってOpenFOAM用メッシュを作成する方法を解説します．snappyHexMeshはSTLファイルからメッシュをヘキサ形状で作成するためのユーティリティです．好きな形状のSTLファイルをCADなどで作成しておき，その形状のメッシュを作成することが可能です．説明上，作業ディレクトリ名を*cylinder*とし，円柱周り流れのメッシュを作成してみます．

動画でも解説してます．
https://www.youtube.com/watch?v=bas-5wbHOoY


完成イメージです．
![screenshot.105.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/551273/36ab329d-5ca7-f56b-4e7d-95acc63d7d47.jpeg)
ブロックの寸法，Lx=0.28m，Ly=0.12m，Lz=0.1mで，原点は円柱中心軸とxy平面の交点としてます．


円柱表面詳細です．境界層も切ります．
![screenshot.106.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/551273/c16d444f-fd74-f308-bc58-1627305d2838.jpeg)

#snappyHexMeshの概要
snappyHexMeshユーティリティはcastllatedMesh，snap，addLayersの3工程で構成されており，処理内容はそれぞれ領域外の除去，境界面の精度向上，境界層用のメッシュ生成，といった感じです．下のようなイメージ．

![screenshot.113.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/551273/1b3f2f34-c298-753b-a567-0015cc065e89.jpeg)


実行には6面体の背景メッシュ，これと矛盾しないstl形式などのジオメトリファイルの2つを用意する必要があります．矛盾している状況の一例として，ジオメトリファイルが背景メッシュより大きい，などが挙げられます．背景メッシュはblockMeshで作成し，ジオメトリファイルはCADなどで作成すればOKです．本記事ではFreeCADで作成します．

この記事はsnappyHexMeshのすべてを解説するものではないので，詳細な仕様を知りたい方は以下を参照ください．
https://cfd.direct/openfoam/user-guide/v7-snappyhexmesh/#x27-1970005.4


**手順**

1. blockMeshで背景メッシュを作成する
2. 円柱のCADファイルを作成する
3. snappyHexMeshでメッシュを作成する

**作業開始時のディレクトリ構成**

以下のようにディレクトリやファイルを用意します．OpenFOAMv7なら
/opt/openfoam7/tutorials/mesh/snappyHexMesh
にあるサンプルディレクトリからひとつ選んでコピペするとcontrolDict，fvSchemes，fvSolution，meshQualityDictを編集しなくて済むので楽です．

```
cylinder/
├── constant/
│   └── triSurface/
└── system/
     ├── blockMeshDict
     ├── controlDict
     ├── fvSchemes
     ├── fvSolution
     ├── meshQualityDict
     └── snappyHexMeshDict
```

#1. blockMeshで背景メッシュを作成する
*cylinder/system/blockMeshDict*を以下のように編集します．

```
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (-0.08 -0.06 0)
    ( 0.2 -0.06 0)
    ( 0.2  0.06 0)
    (-0.08  0.06 0)
    (-0.08 -0.06 0.1)
    ( 0.2 -0.06 0.1)
    ( 0.2  0.06 0.1)
    (-0.08  0.06 0.1)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (100 80 10) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    z_min 
    {
        type            patch;
        faces           ((0 3 2 1));
    }
    z_max
    {
        type            patch;
        faces           ((4 5 6 7));
    }
    y_min
    {
        type            patch;
        faces           ((0 1 5 4));
    }
    y_max
    {
        type            patch;
        faces           ((2 3 7 6));
    }
    x_min
    {
        type            patch;
        faces           ((0 4 7 3));
    }
    x_max
    {
        type            patch;
        faces           ((1 2 6 5));
    }
);

mergePatchPairs
(
);

// ************************************************************************* //
```
パッチ名やメッシュ分割数は特に制限はありません．
*blockMeshDict*を編集し終えたらblockMeshを実行します．

```
cylinder$ blockMesh
```
これで背景メッシュを作成できました．本体は*cylinder/constant/polymesh*にあります．

#2. 円柱のCADファイルを作成する

次にFreeCADで円柱を作成します．わたしはWindowsで作成してます．

直径は0.04mとします．高さはLzを超えてはならないので0.1mとします．ただし，FreeCADからstl形式でエクスポートすると仕様で大きさが1000倍になるので，FreeCAD上で直径0.04mm，高さ0.1mmで作成する必要があります．snappyHexMeshは背景メッシュの座標系とCADファイルの座標系が一致する仕様となっています．円柱の位置をコントロールしたい場合は座標系に注意しながら作成する必要があります．

作成例
![6_完成.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/551273/46351adc-5fa3-ee8a-3fd3-8c92e85b417c.jpeg)

cylinder.stlという名前で出力し*cylinder/constant/triSurface*に置きます．


# 3. snappyHexMeshでメッシュを作成する
*cylinder/system/snappyHexMeshDict*を次のように設定します．

```
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      snappyHexMeshDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

castellatedMesh true;
snap            true;
addLayers       true;

geometry
{
    cylinder
    {
        type triSurfaceMesh;
        file "cylinder.stl";
    }
};

castellatedMeshControls
{
    maxLocalCells 100000;
    maxGlobalCells 2000000;
    minRefinementCells 2;
    nCellsBetweenLevels 4;
    resolveFeatureAngle 30;
    locationInMesh (0.18 0.05 0.09); // Inside point
    allowFreeStandingZoneFaces true;
	
    features
    (
    );

    refinementSurfaces
    {
        cylinder
        {
            level (2 2);
        }
    }

    refinementRegions
    {
    }

}

snapControls
{
    nSmoothPatch 3;
    tolerance 1.0;
    nSolveIter 300;
    nRelaxIter 5;
    nFeatureSnapIter 10;
    implicitFeatureSnap false;
    explicitFeatureSnap true;
    multiRegionFeatureSnap true;
}

addLayersControls
{
    relativeSizes true;
    
	layers
    {
        "cylinder*"
        {
            nSurfaceLayers 5;
        }
    }

    expansionRatio 1.2;
    finalLayerThickness 0.3;
    minThickness 0.25;
    nGrow 0;
    featureAngle 30;
    nRelaxIter 5;
    nSmoothSurfaceNormals 1;
    nSmoothNormals 3;
    nSmoothThickness 10;
    maxFaceThicknessRatio 0.5;
    maxThicknessToMedialRatio 0.3;
    minMedianAxisAngle 90;
    nBufferCellsNoExtrude 0;
    nLayerIter 50;
    nRelaxedIter 20;
}

meshQualityControls
{
    #include "meshQualityDict"

    relaxed
    {
        maxNonOrtho 75;
    }
}



writeFlags
(
    scalarLevels    // write volScalarField with cellLevel for postprocessing
    layerSets       // write cellSets, faceSets of faces in layer
    layerFields     // write volScalarField for layer coverage
);

mergeTolerance 1E-6;

// ************************************************************************* //

```

すべてのパラメータについて説明すると煩雑になるので，経験上重要だと感じたパラメータをかいつまんで上から説明します（というか触ったことないパラメータ結構ある）．

FoamFile直後にあるcastellatedMesh，snap，addLayersは各処理を実行するかどうかを決めるパラメータです．境界層メッシュを切らない場合はaddLayersをfalseにするといった感じで使います．今回は境界層メッシュも切るので全部trueとします．

geometryディクショナリでは読み込むCADファイルを指定します．cylinderはstlファイルの境界（パッチ）の名前のことで，名前は好きにつけてOKです．今回の場合ならパッチとは円柱表面を指します．ここでつけた名前を下の方で使います．typeはtriSurfaceMeshとし，fileにはtriSurfaceに置いたcylinder.stlを指定します．

castellatedMeshControlsは計算領域外の除去処理の設定ディクショナリです．nCellsBetweenLevelsの数字が大きいほど表面周りの解像度が高くなります．わたしは3~4にしてます．locationInMeshはメッシュを作成したい領域のどこかの座標を適当に設定します．つまり背景メッシュの内側かつ円柱の外側のどこかです．円柱の内側を指定すると，円柱内部のメッシュが作成され，外側は除去されてしまいます．refinementSurfacesはパッチごとにを設定するパラメータです．各パッチについて最小，最大(min max)の解像度を指定します．わたしはどちらも2ぐらいにしてます．

snapControlsはsnap処理の設定ディクショナリですが触ったことありません．

addLayersControlはaddLayers処理の設定ディクショナリです．relativeSizesはtrueとすると，境界層の外側の格子に合わせて境界層内の格子幅が計算されます．layersディクショナリではパッチごとの積層数を指定します．パッチ名が"cylinder*"とワイルドカードになってるのは，"[geometryディクショナリで指定したパッチ名]*"としないとちゃんと処理されないからです．私のようなOpenFOAM初学者には分かりにくい部分でした．expansionRatioは層間の格子幅比率です．

編集し終えたらsnappyHexMeshを実行します．

```
cylinder$ snappyHexMesh               // 時間ディレクトリが作成されてその中に出力
cylinder$ snappyHexMesh -overwrite    // polyMeshを上書き出力
```
終わったら名前が数字のディレクトリが新たに3つ作成されます．一番大きな数字のディレクトリ内に作成し終えたメッシュ本体となるpolyMeshディレクトリが入っています．このディレクトリを流体解析用のジョブの*constant*ディレクトリにコピペして使ってください．

終わり．