OpenFOAMで翼形状の非定常計算を行う

## 更新履歴

- 2025/01/22: 初版
- 2025/02/20: 計算環境とOpenFOAM設定を修正
- 2025/04/19: boundaryの設定を修正

## はじめに

本稿では[こちらのチュートリアル](https://www.openfoam.com/documentation/guides/latest/doc/verification-validation-naca0012-airfoil-2d.html#sec-verification-validation-naca0012-airfoil-2d-resources)を参考に、"NACA 0012"翼型モデルの非圧縮性乱流モデルを非定常で計算する備忘を綴る。

## 環境

- AWS Parallel Cluster: ```3.9.20```
- ヘッドノードインスタンスタイプ：```r6g.medium```
- 計算ノードインスタンスタイプ: ```c7gn.xlarge```
- OpenFOAM: ```OpenFOAM-v2306```
- スケジューラ: ```Slurm```

大体これで12hくらい計算時間がかかります。[こちら](https://aws.amazon.com/jp/ec2/pricing/on-demand/)のEC2インスタンスタイプ料金表を基に一回あたりの計算コストを算出：

```math
(0.0608_{[\$]}+0.3148_{[\$]})\times12_{[h]}=4.5072_{[\$]}
```

よって、1計算あたり約700円弱となる。

## チュートリアルケースの取得

以下のチュートリアルケースを作業ディレクトリにコピーしてくる。

```bash
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/airFoil2D .
```

```tree```コマンドでファイルとフォルダを確認すると以下のようになる。

```bash
% tree .
.
├── 0.orig
│   ├── U
│   ├── nuTilda
│   ├── nut
│   └── p
├── Allclean
├── Allrun
├── constant
│   ├── polyMesh.orig
│   │   ├── boundary.gz
│   │   ├── cells.gz
│   │   ├── faces.gz
│   │   ├── neighbour.gz
│   │   ├── owner.gz
│   │   └── points.gz
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── decomposeParDict
    ├── fvSchemes
    └── fvSolution

4 directories, 18 files
```

主な各ファイルの役割は次の通り。

1. ```0.orig/``` 
    場の変数の設定ファイル（速度U、圧力pなど）が含まれるディレクトリ
1. ```Allrun```
    1. 境界条件や初期値の設定を行う```0```を```0.orig```から、そして```constant/polymesh```を```constant/polymesh.orig```から参照
    1. 並列計算を実行
1. ```constant/```
    メッシュや物性値、乱流モデルなどの設定が含まれるディレクトリ
    1. ```polyMesh/```
        メッシュ情報を含むディレクトリ。チュートリアル内でplot3dToFoamを用いたメッシュ分割結果が格納
    1. ```transportProperties/```, ```turbulenceProperties/```
        物性値の設定ファイル
1. ```system/```
    計算の制御の設定ファイルなどが含まれるディレクトリ
    1. ```controlDict/```
        計算の制御の設定ファイル。計算の終了時刻や時間刻み幅、結果の出力タイミングなどを設定
    1. ```fvSolution/```
        代数方程式ソルバーの設定と、SIMPLE法などの設定を含むファイル
    1. ```fvSchemes/```
        離散化スキームの設定ファイル

### Paraviewによるメッシュの確認

作業ディレクトリ直下に```case.foam```などを作成して、ParaViewでメッシュを確認する。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/36017025-8567-d4d5-6c80-6a475595ea35.png)

## メッシュの置き換え

メッシュは上のOpenFOAMチュートリアルで用意されているものではなく、NASAが提供するNACA 0012翼型に特化した計算格子があるので、それを用いる。メッシュの詳細については[こちら](https://turbmodels.larc.nasa.gov/naca0012_grids.html)を参考

本稿でもそれに倣って、[(2-D) 225 x 65 (129 points on airfoil surface)](https://turbmodels.larc.nasa.gov/NACA0012_grids/n0012_225-65.p2dfmt)を用いる。作業ディレクトリ以下で以下のコマンドを実行

```bash
cd constant/
wget https://turbmodels.larc.nasa.gov/NACA0012_grids/n0012_225-65.p2dfmt
cd ..
plot3dToFoam constant/n0012_225-65.p2dfmt -2D 1 -noBlank
```

```constant/polyMesh/```が生成されているので、改めてメッシュを可視化して確認してみる。

**【遠距離ビュー】**
![far.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/74da517c-eb1c-c046-3b36-6117ab219622.png)

**【近距離ビュー】**
![near.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/ad0a1221-79b3-5d17-fc8d-f6c2e6c5a787.png)

### パッチの分割

上記の操作により作成された```constant/polyMesh/boundary```を見ると```defaultFaces```のみになっているので、以下のコマンドでパッチを自動分割する。

```shell
autoPatch -overwrite 30
```

こうすると```auto_*```というパッチ名が自動で割り当てられているので、分かり易く次のように書き換える。面倒であるが、ParaViewを使って1パッチずつ対応を確認した。

```cpp: constant/polyMesh/boundary
5
(
    back
    {
        type            empty;
        physicalType    empty;
        nFaces          14336;
        startFace       28432;
    }

    front
    {
        type            empty;
        physicalType    empty;
        nFaces          14336;
        startFace       42768;
    }

    outlet
    {
        type            patch;
        physicalType    outlet;
        nFaces          128;
        startFace       57104;
    }

    aerofoil
    {
        type            wall;
        physicalType    wall;
        nFaces          128;
        startFace       57232;
    }
    
    inlet
    {
        type            patch;
        physicalType    inlet;
        nFaces          224;
        startFace       57360;
    }
)
```

パッチが定義できたところで以下のコマンドで```polyMesh.orig```を置き換える。

```bash
rm -rf constant/polyMesh.orig/*
cp -r constant/polyMesh/* constant/polyMesh.orig/
```

## 空力係数の計算

```forceCoeffs```関数を用いて```Cl```値/```Cd```値など空力係数を計算する。

```cpp: system/forceCoeffs
forces
{
    type            forceCoeffs;
    libs            ( forces );
    writeControl    timeStep;
    writeInterval   1;
    patches         ( aerofoil );
    rho             rhoInf;
    log             true;
    rhoInf          1;
    liftDir         ( 0 1 0 );
    dragDir         ( 1 0 0 );
    CofR            ( 0.25 0 0 ); // 回転中心座標。翼先頭から1/4地点
    pitchAxis       ( 0 0 1 );
    magUInf         51.4815;
    lRef            1; // 翼の長さ (NASA 2DグリッドでX方向=1)
    Aref            1; // =翼の長さ×スパン長（NASA 2DグリッドでZ方向=1）
}
```

## PISO法による非定常計算

OpenFOAMには非定常解析ソルバーとして、PISO法やPIMPLE法などが用意されている。本章ではPISO法アルゴリズム概要について簡単に記載して置く。

### ナビエ・ストークス（Navier-Stokes/NS）方程式

一般に、流速がマッハ数0.3以下（密度変化が5%以下）であれば非圧縮性流れとして考える。後節で記述するが、今回の翼形状周りの流れシミュレーションには流速$=51.4815_{[m/sec]}\simeq0.15_{Ma}$を与えるので、非圧縮性方程式を考える。

流速を$\boldsymbol{u}$、圧力を$\boldsymbol{p}$、流体密度を$\rho$、粘性応力テンソルを$\boldsymbol{\tau}$として、ナビエ-ストークス方程式は次式で表される：

```math
\cfrac{\partial\boldsymbol{u}}{\partial t}
+
\left(
    \boldsymbol{u}\cdot\nabla
\right)
\boldsymbol{u}
=
-\cfrac{1}{\rho}\nabla\boldsymbol{p}
+
\nabla\cdot\boldsymbol{\tau}
```

ここで、粘性応力テンソル$\boldsymbol{\tau}$はせん断粘性係数を$\mu$、粘性係数を$\lambda$として以下のように分解される：

```math
\boldsymbol{\tau}
=\mu
\lbrack
    \nabla\boldsymbol{u}
    +(\nabla\boldsymbol{u})^\mathsf{T}
\rbrack
+\lambda
(
    \nabla\cdot\boldsymbol{u}
)
\boldsymbol{1}
```

さらに多くの場合、体積粘性率$\chi=\lambda+(2/3)\mu$はストークス近似のもとで0とされる。
故に動粘性係数を$\nu:=\mu/\rho$として元の式に代入すると：

```math
\cfrac{\partial\boldsymbol{u}}{\partial t}
+
\left(
    \boldsymbol{u}\cdot\nabla
\right)
\boldsymbol{u}
=
-\cfrac{1}{\rho}\nabla\boldsymbol{p}
+
\nabla\cdot
\lbrack
    \nu\lbrace
        \nabla\boldsymbol{u}
        +(\nabla\boldsymbol{u})^\mathsf{T}
    \rbrace
\rbrack
-
\nabla
\left(
    \frac{2}{3}\nu\nabla\boldsymbol{u}
\right)
```

### PISO法による圧力-流速連成

OpenFOAMのPISO法では以下の計算手順により圧力と流速を連成させる：

1. **予測速度場の計算**
   ナビエ-ストークス方程式を空間離散化して次のような式を構築する：

   ```math
   A_P\boldsymbol{u}_P + \sum A_N\boldsymbol{u}_N = -\nabla p
   ```

   ここで、$A_P$および$A_N$はそれぞれ注目セルと隣接セルの係数である。
   これを次の形に分離する：

   ```math
   A\boldsymbol{u} = \boldsymbol{H} - \nabla p
   ```

   圧力項を一旦無視して予測速度場$\hat{\boldsymbol{u}}$を計算する：

   ```math
   A\hat{\boldsymbol{u}} = \boldsymbol{H}
   ```

   OpenFOAMにおける`momentumPredictor`オプションは、この仮速度場の計算を実行するスイッチである。

2. **圧力修正方程式の構築**
   非圧縮性条件$\nabla\cdot\boldsymbol{u}=0$を満たすため、圧力補正方程式を構築する。
   仮速度場$\hat{\boldsymbol{u}}$を基に、次のような式が導かれる：
   
   ```math
   \nabla\cdot\left(\frac{1}{\rho}\nabla\boldsymbol{p}'\right) = \nabla\cdot\hat{\boldsymbol{u}}
   ```

   ここで、$\boldsymbol{p}'$は圧力補正項を表す。

3. **圧力場の更新**
   圧力補正項$\boldsymbol{p}'$を用いて圧力場を更新する：
   
   ```math
   \boldsymbol{p}^{(n+1)} = \boldsymbol{p}^{(n)} + \boldsymbol{p}'
   ```

4. **速度場の修正**
   修正された圧力場$\boldsymbol{p}'$を用いて速度場を次のように修正する：

   ```math
   A\boldsymbol{u}^{(n+1)} = \boldsymbol{H} - \nabla\boldsymbol{p}^{(n+1)}
   ```

   この修正はPISO法の反復処理中に`nCorrectors`オプションによって指定された回数だけ行われる。

5. **残差の監視**
   各計算ステップの残差（例えば速度、圧力の収束状況）は`solverInfo.dat`ファイルに出力される。このファイルには、各反復ステップの残差値が記録され、収束基準を確認するために使用する。

6. **反復処理**
   圧力補正と速度修正を複数回繰り返すことで、連続の式$\nabla\cdot\boldsymbol{u}=0$をより高い精度で満たすように調整する。必要に応じて、圧力補正方程式を再構築し、精度を向上させる。

## 非定常計算への変更

コピーした定常計算チュートリアルとNASA提供のグリッドを元に非定常計算を行う変更を行う。

### ```controlDict```

前述の通り、今回は代数方程式ソルバーに```pisoFoam```を用いる。

```cpp: system/controlDict
application     pisoFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         5.0;
deltaT          1e-4;
writeControl    runTime;
writeInterval   1e-4;
purgeWrite      2001;
writeFormat     ascii;
writePrecision  10;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;
functions
{
    #include "forceCoeffs"
    #includeFunc solverInfo(U,p,k,omega)   
}
```

#### ```deltaT```と```writePrecision```について

非定常シミュレーション結果として、システムの振動が観測できる値にする。例えば、以下のようにグラフ可視化しても、データを直接表示しても同じ値を取るケースなどが生じうる。

**【```postProcessing/forces/0/coefficient.dat```の描画結果】**

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/1569b3db-f6ed-aff0-dd5a-118345271d8a.png)

**【```postProcessing/forces/0/coefficient.dat```のデータ中身】**

![image (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/2d434acc-b604-10ae-469b-110afc2e1c91.png)

このように、系の振動を観測できない場合は```deltaT```を小さく、```writePrecision```を大きくするのを試してみると良さそうである。

#### ```endTime```について

後述する仰角```alphaDeg```にもよるが、```coefficient.dat```を観察してみると大5秒もあれば収束することが分かる。あるいは、残差```solverInfo.dat```の収束を見て設定しても良い。
    
### ```fvScheme```

以下の点だけ修正を加える。```backward```で後方二次精度となる。

```cpp
ddtSchemes
{
    default         backward;
}
```

### ```fvSolution```

| Field   | Linear Solver     | Smoother              | Tolerance (rel) |
|---------|-------------------|-----------------------|-----------------|
| U       | Smooth solvers    | Gauss Seidel Smoother | 0.01            |
| p       | GAMG Solver       | Gauss Seidel Smoother | 0.01            |
| nuTilda | Smooth solvers    | Gauss Seidel Smoother | 0.01            |

#### 数値スキーム設定```system/fvSolution```

```cpp: system/fvSolution
solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-09;
        relTol          1e-02;
    }
    pFinal
    {
        $p;
        relTol          0;
    }
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        relTol          1e-02;
        tolerance       1e-09;
    }
    nuTilda
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        relTol          1e-02;
        tolerance       1e-09;
    }
}

PISO
{
    momentumPredictor        yes;
    nCorrectors              3;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        nuTilda         0.7;
    }
}
```

#### ```tolerance```について

計算結果として出力される```postProcessing/solverInfo(U,p,k,omega)/0/solverInfo.dat```を見てみると、以下のようなグラフを可視化することができる。

**【```postProcessing/solverInfo(U,p,k,omega)/0/solverInfo.dat```の描画結果】**

![angle_0_residuals.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/0317fdd2-121f-685a-575d-b104a5288ca4.png)

まだ十分に収束しているとは言えないが、設定した$1e-8$付近のオーダーで残差が推移していることが分かる。
加えて、この時の```coefficient.dat```を併せて見てみると

**【```postProcessing/forces/0/coefficient.dat```の描画結果】**

![angle_0_coeffs.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/dab0ee31-17a8-e091-ec7b-d4216b39793c.png)

このように観測データである空力係数```Cl```, ```Cd```が一定時間経過後に収束していることが読み取れる。（前述した振動はグラフでは見えないレベルで生じている）一方で、```torelance```の値が大き過ぎる、或いは設定しないと以下のように空力係数の値が収束（厳密には微振動）せず、大きな振れ幅のまま発散してしまう。

**【空力係数が発散するケース】**

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/255c8eb5-1918-32dd-369f-346db5151c0d.png)

実際、このようなケースは後述する```Cd```, ```Cl```プロットで比較データと大きく差異が生じてしまったため、良いとは推測し難い。このように、空力係数と残差の時間推移を観察しながら適切な```torelance```の決定することで、PISO法による代数方程式計算の打ち切り（精度）を調整する。

### 乱流モデル```constant/turbulenceProperties```の設定

2次元シミュレーションでのLESモデルの使用に関して、LESモデルは通常、3次元非定常計算で用いられるが、2次元問題や定常計算でも選択可能となっている。ところが、2次元でのLESモデルの適用は物理的に適切でないため、今回はSpalart-Allmarasモデルを採用する。参考は[こちら](https://www.cfd-online.com/Forums/openfoam/121511-les-2d-case.html)。

```cpp:constant/turbulenceProperties
simulationType          RAS;

RAS
{
    RASModel            SpalartAllmaras;
    turbulence          on;
    printCoeffs         on;
}
```

## 場の設定

チュートリアルを再現するために、場の設定を変更する。

### ```Allrun```

今後の拡張として仰角```alphaDeg```を変えてCFDを流したいので、```foamDictionary```コマンドを利用して場や制御の変数を書き換えるようにする。

```bash: Allrun
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions

#------------------------------------------------------------------------------
# Variables

parallel=true # flag to enable computations in parallel mode
alphaDeg=0 # Angle in degrees (adjust as needed)
speed=51.4815 # Speed (adjust as needed)

# Calculate cos(α) and sin(α)
alphaRad=$(echo "$alphaDeg * 3.141592653589793238462643 / 180" | bc -l)
cosAlpha=$(echo "c($alphaRad)" | bc -l)
sinAlpha=$(echo "s($alphaRad)" | bc -l)

#------------------------------------------------------------------------------
# Set initial speed field

# Modify `0.orig/U` internalField
foamDictionary 0.orig/U \
    -entry internalField \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

# Modify `0.orig/U` boundaryField inlet
foamDictionary 0.orig/U \
    -entry boundaryField.inlet.freestreamValue \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

# Modify `0.orig/U` boundaryField outlet
foamDictionary 0.orig/U \
    -entry boundaryField.outlet.freestreamValue \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

#------------------------------------------------------------------------------
# Direction setting of force coefficient (lift and drag)

# Modify `system/forceCoeffs` liftDir
foamDictionary system/forceCoeffs \
    -entry forces.liftDir \
    -set "($(
        echo "-1 * $sinAlpha" | bc -l
    ) $(
        echo "$cosAlpha" | bc -l
    ) 0)"

# Modify `system/forceCoeffs` dragDir
foamDictionary system/forceCoeffs \
    -entry forces.dragDir \
    -set "($(
        echo "$cosAlpha" | bc -l
    ) $(
        echo "$sinAlpha" | bc -l
    ) 0)"

#------------------------------------------------------------------------------
# Initialization

restore0Dir
[ -d constant/polyMesh ] || cp -rf constant/polyMesh.orig constant/polyMesh

#------------------------------------------------------------------------------
#  Running Parallel Computations

if [ "$parallel" = true ]
then
    runApplication decomposePar
    runParallel $(getApplication)
else
    runApplication $(getApplication)
fi
```

### ```0.orig```

#### 速度場```U```

| Patch        | Condition          | Value [$m\cdot s^{-1}$]       |
|--------------|--------------------|---------------------|
| Inlet        | freestreamVelocity | 51.4815                 |
| Outlet       | freestreamVelocity | 51.4815                 |
| Sides (y-dir)| empty              | -                  |
| Aerofoil     | fixedValue         | (0.0, 0.0, 0.0)    |

上記の```Allrun```により```internalField```フィールドは書き換えられるが、安全のために```51.4815```で初期化しておく。

```cpp: 0.orig/U
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (51.4815 0 0);

boundaryField
{
    inlet
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}
```

#### 圧力場```p```

| Patch        | Condition          | Value [$m^2\cdot s^{-2}$]       |
|--------------|--------------------|---------------------|
| Inlet        | freestreamPressure | 0    |
| Outlet       | freestreamPressure | 0    |
| Sides (y-dir)| empty              | -    |
| Aerofoil     | zeroGradient       | -    |

```cpp: 0.orig/p
dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            freestreamPressure;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestreamPressure;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            zeroGradient;
    }
}
```

#### 乱流動粘性係数```nut```

| Patch        | Condition          | Value [$m^2\cdot s^{-1}$]       |
|--------------|--------------------|---------------------|
| Inlet        | freestream | $8.58\times10^{-6}$    |
| Outlet       | freestream | $8.58\times10^{-6}$    |
| Sides (y-dir)| empty              | -    |
| Aerofoil     | fixedValue       | 0    |

```cpp: 0.orig/nut
dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 8.56e-6;

boundaryField
{
    inlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform 0.0;
    }
}
```

#### 乱流動粘性```nuTilda```

| Patch        | Condition          | Value [$m^2\cdot s^{-1}$]       |
|--------------|--------------------|---------------------|
| Inlet        | freestream | $3.432\times10^{-5}$    |
| Outlet       | freestream | $3.432\times10^{-5}$    |
| Sides (y-dir)| empty              | -    |
| Aerofoil     | fixedValue       | 0    |

```cpp: 0.orig/nuTilda
dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 3.432e-5;

boundaryField
{
    inlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform 0.0;
    }
}
```

## シミュレーション結果

本シミュレーション結果の検証には、NACA 0012翼型の揚力係数```Cl```および抗力係数```Cd```と迎角（Angle of Attack）```alphaDeg```との関係を示したデータセットを使用した。このデータセットは、以下のような実験的/数値的な結果に基づいく。

- 実験データ: [Ladsonらによる実験データ](https://turbmodels.larc.nasa.gov/NACA0012_validation/CLCD_Ladson_expdata.dat)
    
    Ladsonらの実験データは、NASA Langley 研究センターで取得されたもので、風洞試験による高精度な```Cl```と```Cd```の測定値を含む。

- 数値データ: [数値シミュレーション (CFL3D ソフトウェア)](https://turbmodels.larc.nasa.gov/NACA0012_validation/n0012clcd_cfl3d_sa.dat)
    NASAが開発したCFL3Dソフトウェアを使用して計算された数値データ。この数値解は、実験データとの整合性を評価するための基準として使用される。

本シミュレーションでは、実験データとの比較を行い、揚力係数および抗力係数の計算結果がこれらの既知データとどの程度一致しているかを検証した。特に、迎角の増加に伴う非線形挙動や失速時の特性を検証の焦点に置く。

### 1. 揚力係数：```Cl```値

![visualize_cl_summary.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/53eb0f79-c07d-b0ac-d339-0ebd0f267f05.png)

### 2. 抗力係数：```Cd```値

![visualize_cd_summary.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/b2bd9f93-76ba-0895-72da-686e9779546d.png)

このように揚力係数```Cl```値、抗力係数```Cd```値いずれのケースにおいても実験データの特徴を良く再現出来ているように見える。仰角```alphaDeg```=16付近を頭打ちにそれ以上の仰角だと、翼状後流の剥離が顕著となり、揚力/抗力共にドラスティックに変化する。

## 流れ場の可視化

ParaViewを用いて翼状まわりの流れ場を可視化する。手順は以下の通り：

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/6bd99ce0-57e9-c4a7-ca70-4b129b32a219.png)

- ```processor*/```ディレクトリがある同じ階層に```case.foam```を作成し、ParaViewでそれを読み込む。
- 画像内①にあるCase Typeを"Decomposed Case"に設定
- 画像内②にあるScalar Sizeを64bitに設定
- Applyを押下
- 画像内④で描画する場を選択

### 仰角```alphaDeg```ごとの流速場分布

#### 1. ```alphaDeg```=0

![angle_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/94be2aa2-82a6-7440-4aa9-83d864cdd32d.png)

#### 2. ```alphaDeg```=10

![angle_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/bbb3e1f5-3499-b371-26ca-4eb6258af7bc.png)

#### 3. ```alphaDeg```=15

![angle_15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/1b3971f2-b859-219e-e61d-199bc9839c53.png)

### 剥離現象の観察

ParaViewの視覚化手法を"Surface LIC (Line Integral Convolution)"を選択して、場の流れを直感的に視覚化してみる。```alphaDeg```=19と仰角が大き過ぎて揚力が落ち込んでいるケースで試してみる。

![paraview.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/409783/0145abe3-7271-9f11-2fa6-9cfaf45ebf17.png)

## おわりに

いい感じに再現できた。OpenFOAMは情報共有がネット上にもあまりなく、苦戦したが本稿を作成するに際して一番役に立ったのは[この30秒のYouTube動画](https://youtu.be/9XKWyihRYnY?si=8CFRyPgBt_TxSMkD)だった。意外とYouTubeで調べて見るのもアリだなと思った。