OpenFOAM+cfMeshを使用してモデルロケット周りの流体解析(1)

#用いたケースファイル
https://github.com/reuscaelum/CFD_cases
これを解説していきます。

#実行環境
- ubuntu 18.04
- OpenFOAM 6
- cfMesh-Development

#用いたCADファイル
![Screenshot from 2019-09-12 23-55-57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/300702/39d26ba5-0c0e-650d-06b6-031c04a098c5.png)
今回の記事用にモデルロケット風なCADデータを用意しました.
機体の進行方向が+Z軸方向です.

#解析領域の作成
cfMeshのsurfaceGenerateBoundingBoxコマンドを使用します.(例:surfaceGenerateBoundingBox rocket.stl rocket_box.stl 1.5 1.5 4 4 20 20)
入力するSTLファイルと出力するSTLファイルを指定して,その後に解析領域の-X方向,+X方向,-Y方向,+Y方向,-Z方向,+Z方向(単位はm)の大きさを指定します.物体の後流の流れに影響を及ぼさない程度に大きくしておいたほうがいいらしいです.今回は過剰に取りすぎた感があります.

#メッシュ生成
今回使用したmeshdict

```:system/meshdict
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                |
| \\      /  F ield         | cfMesh: A library for mesh generation          | 
|  \\    /   O peration     |                                                |
|   \\  /    A nd           | Author: Franjo Juretic                         | 
|    \\/     M anipulation  | E-mail: franjo.juretic@c-fields.com            |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version   2.0;
    format    ascii;
    class     dictionary;
    location  "system";
    object    meshDict;
}

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

surfaceFile "rocket_box.stl";

maxCellSize 5.0;

boundaryCellSize 0.1;

minCellSize 3;



localRefinement
{
    "modelrocket.*"
    {
        cellSize 0.0005;
    } 
}

boundaryLayers
{
  //nLayers 15;

  //thicknessRatio 1.2;

  //maxFirstLayerThickness 0.5;

    patchBoundaryLayers
    {
        "modelrocket.*"
        {
            nLayers           10;

            thicknessRatio    1.02;

            maxFirstLayerThickness 0.0001;

            allowDiscontinuity 0;
        }



    }
}

renameBoundary
{
    defaultName		fixedWalls;
    defaultType		wall;

    newPatchNames
    {
        "modelrocket.*"
        {
            newName	rocket;
	    type     wall;
        }
 
        "zMax"
        {
            newName     inlet;
            type     patch;
        }

        "zMin"
        {
            newName     outlet;
            type     patch;
        }
    }
}
               
     
// ************************************************************************* //
```

systemの中にあるmeshdictでメッシュの設定をします.まずsurfaceFileで読み込むstlファイルを指定します.この指定するstlは先程解析領域で出力したstlファイルを指定します.
そのあとにメッシュのサイズなどを設定する項目がありますが飛ばします.

localRefinementの項目でstlファイル周辺のメッシュを細分化します.今回は周辺のメッシュを0.5mmに細分化します.
次にboundaryLayersの項目で境界層の設定をします.CFD OnlineのY+ Estimate https://www.cfd-online.com/Tools/yplus.php
などを利用して十分な境界層メッシュが生成されるようにします.今回は乱流モデルとしてK-ω SSTを用いるので10層の境界層メッシュを挿入します.
次に境界層メッシュの成長率と1層目の最大厚さ,不連続なメッシュを許容するかという設定をします.
次に境界の名前と性質を変更します.流入面,流出面を定義しておきます.
meshdictを保存したらcartesianMeshを実行してmeshを生成します. 
checkMeshをしてメッシュのアスペクト比や歪度などを見てダメそう(アスペクト比が200超とか)だと設定を変えてメッシュをもう一度生成したほうが良さそうです.
またparaViewで狙い通りのメッシュが生成されているかも見ておいたほうが良さげです.数字上では良くても,壊れているメッシュが生成されていたりします.(今回は大丈夫でした)

今回生成したメッシュ
![Qiita_mesh.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/300702/047eb6d9-4f52-6106-b4ca-7e375b9dd817.png)

#乱流モデルの選択&初期値の設定
乱流モデルはconstantの中のturbulence Propertiesで選択します.今回はk-ω SSTを用います.詳しい説明はこちらをどうぞhttps://turbmodels.larc.nasa.gov/sst.html

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

simulationType RAS;

RAS
{
    RASModel            kOmegaSST;

    turbulence          on;

    printCoeffs         on;
}

// ************************************************************************* //
```
次に初期値を設定します.今回は-Z方向に100[m/s]の空気が流入,流出面で大気圧開放するようにします.
tutorialのmotorbikeのファイルを持ってきて今回のケースに合うように境界の名前を書き換えたりします.
今回はK-ω SSTを用いるのでUとpの他にk,nut,omegaの初期値を設定します.

```:0/U
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      binary;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    rocket
    {
        type            noSlip;
    }
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 -100);
    }
    outlet
    {
        type            zeroGradient;
    }
    fixedWalls
    {
        type            noSlip;
    }
}


// ************************************************************************* //
```
```:0/p
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      binary;
    class       volScalarField;
    location    "0";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    rocket
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    fixedWalls
    {
        type            zeroGradient;
    }
}


// ************************************************************************* //
```

次にsystemのcontroldictで解析の設定をします.
applicationでどのソルバーで解析するかを設定します.今回はsimpleFoamを選択します.
次に何回反復計算するかを設定します.simpleFoamを用いるのでendTimeだけ気にすればいいと思います.今回は1000回反復計算させてみます.deltaTは今回は定常計算なので1にします.
次にfunctionsで物体にかかる力のログ出力の設定をします.
またforceCoeffsを読み込んで抗力係数を算出します.

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

application     simpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         1000;

deltaT          1;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     binary;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


functions
{
forces
{  
type	forces; 
functionObjectLibs ("libforces.so"); 

outputControl	timeStep;
timeInterval 	1;
log	yes;
patches	(rocket);
p	p;
U	U;
rho	rhoInf;
log	true;
rhoInf  1;   
CofR	(0 0 0);
}
    #include "forceCoeffs"
}

// ************************************************************************* //
```
次に離散化スキームと代数方程式ソルバーの設定をfvSchemesとfvSolutionで設定します.ここは自分が説明するよりPENGUINITSさんのページを見るのが早いと思うので見てくださいhttp://penguinitis.g1.xrea.com/study/OpenFOAM/index.html

#実行させるシェルスクリプトを準備
tutorialからAllrunとAllcleanをコピペして使います.decomposeParを使って並列計算させています.メッシュをみたいだけ用にAllmeshというシェルスクリプトも用意しました.

```:Allrun
#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

surfaceGenerateBoundingBox rocket.stl rocket_box.stl 1.5 1.5 4 4 20 20
runApplication cartesianMesh
runApplication checkMesh -constant
runApplication decomposePar
runParallel $(getApplication)

runApplication reconstructParMesh -constant
runApplication reconstructPar -latestTime
runApplication foamLog log.simpleFoam
#------------------------------------------------------------------------------
```
```:Allclean
#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial clean functions
. $WM_PROJECT_DIR/bin/tools/CleanFunctions

# Remove surface and features


cleanCase

#------------------------------------------------------------------------------
```

加えてこの上の階層に次のようなシェルスクリプトを置いて複数のケースを一回のシェルスクリプトの実行で回せるようにします.

```:start
#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

dir_path="${0%/*}"
dirs=`find $dir_path -maxdepth 1 -type d -not -name .`;

for dir in $dirs;
do
    
    cd $dir
    ./Allrun
    cd ..
done
#------------------------------------------------------------------------------
```
これを実行すると解析が始まります.

#結果
ノーズ後端で剥離してるのがなんとなくわかりますね
![result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/300702/0c49a4fa-d5aa-7969-4464-010487e41abf.png)

#感想
技術系記事を書くのが初めてなので至らないところがあったら申し訳ないです.全国のロケッティアが参考になればと思い書いてみました.圧縮性のソルバーもrunApplicationの指定変えて,pの初期値の次元を合わせたらできるので超音速機作る際はそちらを使えばいいと思います.
質問があればtwitter:@8900_blueskyまでどうぞ
次はスペースバルーンの軌道シミュレーションの記事についてor圧縮性ソルバーを用いた解析について書こうと思います.(もっとスペースバルーンやるサークル増えて…)




#参考にしたページ
https://damassets.autodesk.net/content/dam/autodesk/www/campaigns/sim-day-2015/Autodesk_CFD_%E8%A7%A3%E6%9E%90%E3%83%86%E3%82%AF%E3%83%8B%E3%83%83%E3%82%AF.pdf
http://penguinitis.g1.xrea.com/study/OpenFOAM/index.html
https://turbmodels.larc.nasa.gov/sst.html