Notched Disc In Rotating Flow

# はじめに
ESI版OpenFOAM v2406を対象にして、タイトルに記した計算を行う。
この計算は、回転する流速場にあるノッチ付き円盤状の水(以下、$\alpha$と記す)が一周してどのくらい初期の形状を保てるか、というものである。
ソルバはinterIsoFoam/interFoamを使用し、$\alpha$の移流のみ計算する。
ESI版はinterIsoFoamのチュートリアルに本問題があり、それを使用する。
Foundation版OpenFOAM 11でも計算を行う。

# 解析モデル
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/34a2439c-fd44-6341-2ad4-956726c31662.png)

<table>
  <caption>設定値一覧</caption>
  <thead>
    <tr>
      <th>名前</th> <th>設定値</th>
    </tr>
  </thead>
  <tr>
    <td> ドメイン長さ $L [m]$ </td> <td>$1$</td>
  </tr>
  <tr>
    <td> ディスク半径 $R [m]$ </td> <td>$0.15$</td>
  </tr>
  <tr>
    <td> ノッチ幅 $W [m]$ </td> <td>$0.05$</td>
  </tr>
  <tr>
    <td> ノッチ高さ $H [m]$ </td> <td>$0.25$</td>
  </tr>
  <tr>
    <td> ディスク初期位置 $(x_0,y_0,z_0)$ </td> <td>$(0.5, 0, 0.75)$</td>
  </tr>
</table>

# 回転する流速場
ESI版OpenFOAMでは、function objectで流速場の初期設定ができる。
![U0_v2406.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/9654a05a-e213-01ca-80ec-b8e636a8a06e.png)

速度場$U(u,0,v)$の分布は次式の通りである。
```math
u = \frac{\pi}{3.14} (0.5 - z)
```
```math
v = \frac{\pi}{3.14} (x - 0.5)
```
```controlDict
functions
{
    setVelocity
    {
        type            setFlow;
        libs            (fieldFunctionObjects);
        writeControl    writeTime;
        mode            rotation;
        scale           1;
        //reverseTime     1;
        omega           1.000507215; // 2*pi/(2*3.14)
        //omega           6.28318530718;
        origin          (0.5 0 0.5);
        refDir          (1 0 0);
        axis            (0 1 0);
    }
}
```
system/controlDictのオリジナルからreverseTimeをコメントにし、omegaを$\frac{2\pi}{2✕3.14}$に変更する。つまり、ノッチ付き円盤の周期は、$2✕3.14=6.28 [s]$である。これを計算終了時刻とする。

# ノッチ付き円盤の動き
![frame.0000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/4f50b03f-c4ef-a732-a61a-52faa60ff581.png)
![frame.0002.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/1eea9c09-99b3-253d-d021-93ef6e8397e8.png)
![frame.0004.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/8925f5c3-9298-964d-2a02-e143af8c850a.png)
![frame.0006.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/6f370fae-61f8-c765-52d9-84382c12d665.png)
![frame.0008.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/2969e4b0-989d-64ea-990e-09d7e6ed7e2f.png)

# 評価式
計算開始時と計算終了時の$\alpha$の分布より、次の誤差を計算する。
```math
\epsilon_{ND2} = \frac{1}{L} \sum_{j=1}^{N} A_{j} |\alpha_{J}^{I} - \alpha_{j}^{F}|
```
<table>
  <caption>変数一覧</caption>
  <thead>
    <tr>
      <th>変数名</th> <th>説明</th>
    </tr>
  </thead>
  <tr>
    <td>  $L$ </td> <td>界面の長さ $[m]$</td>
  </tr>
  <tr>
    <td> $A_j$ </td> <td>$セルjの面積 [m^2]$</td>
  </tr>
  <tr>
    <td> $\alpha_j^I$ </td> <td>$計算開始時のセルjの\alpha$</td>
  </tr>
  <tr>
    <td> $\alpha_j^F$ </td> <td>$計算終了時のセルjの\alpha$</td>
  </tr>
</table>

$L$は計算開始時の界面の長さとすると、約$1.442 [m]$である。$A_j$はセルの分割数から算出した値とする。文献では$\alpha_j^I$を1タイムステップ実行後の値としているが、ここでは初期値とする。

# ESI版OpenFOAMの計算結果
<table>
  <caption>誤差$\epsilon_{ND2}$の計算結果</caption>
  <thead>
    <tr>
      <th></th> <th>ケース1</th> <th>ケース2</th> <th>ケース3</th>
    </tr>
  </thead>
  <tr>
    <td> 分割数 </td> <td>$100✕100$</td> <td>$200✕200$</td> <td>$400✕400$</td>
  </tr>
  <tr>
    <td> $interIsoFoam$ </td> <td>$0.001 97$</td> <td>$0.000 51$</td> <td>$0.000 24$</td>
  </tr>
  <tr>
    <td> $interFoam$ </td> <td>$0.003 42$</td> <td>$0.001 39$</td> <td>$0.000 75$</td>
  </tr>
</table>

なお、interIsoFoamのreconstructSchemeはisoAlphaである。

# Foundation版の回転する速度場
Foundation版では、0/UのinternalFieldsにcodeStreamを適用することで、初期設定ができる。

```0/U
internalField   #codeStream
{
    codeInclude
    #{
        //#include "fvCFD.H"
        #include "volFields.H"
    #};

    codeOptions
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};

    codeLibs
    #{
        -lfiniteVolume \
        -lmeshTools
    #};

    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());

        vectorField U_0(mesh.nCells(), Foam::Vector<double>(0.0,0.0,0.0));

        forAll(U_0, i)
        {
            const scalar x = mesh.C()[i][0];
            const scalar z = mesh.C()[i][2];
            const scalar u = constant::mathematical::pi * (0.5 - z) / 3.14;
            const scalar v = constant::mathematical::pi * (x - 0.5) / 3.14;
            U_0[i] = Foam::Vector<double>(u,0.0,v);
        }

        // ESI style U_0.writeEntry("", os);
        writeEntry(os, "", U_0);
    #};
};
```
バージョン11以前では、codeIncludeで"fvCFD.H"をインクルードする。
また、ESI版で使う場合は、writeEntryの書き方を変える。
このコードは実行時にコンパイルされ、速度場は次のようになる。
![U0_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/28b3b5eb-575f-5240-3f00-98631bedfe56.png)

この速度場は、ESI版と同じ（のはず）である。

# Foundation版OpenFOAMの計算結果
Foundation版のinterFoamでは、$\alpha$の発散スキームとして次の3つが用意されている。

<table>
  <caption>誤差$\epsilon_{ND2}$の計算結果</caption>
  <thead>
    <tr>
      <th></th> <th>ケース1</th> <th>ケース2</th> <th>ケース3</th>
    </tr>
  </thead>
  <tr>
    <td> 分割数 </td> <td>$100✕100$</td> <td>$200✕200$</td> <td>$400✕400$</td>
  </tr>
  <tr>
    <td> $interfaceCompression$ </td> <td>$0.00411$</td> <td>$0.00177$</td> <td>$0.00097$</td>
  </tr>
  <tr>
    <td> $PLIC$ </td> <td>$0.00449$</td> <td>$0.00177$</td> <td>$0.00097$</td>
  </tr>
  <tr>
    <td> $MPLIC$ </td> <td>$0.00449$</td> <td>$0.00210$</td> <td>$0.00133$</td>
  </tr>
</table>

# まとめ
本問題に関しては、ESI版の方が誤差が小さく、計算の準備が容易である。