OpenFOAMにおける乱流境界条件の設定について

# はじめに

乱流パラメーター(乱流エネルギーや消散速度など)や境界条件の設定について解説します。
すでにPenguinitisさんの[境界条件の設定のページ](http://penguinitis.g1.xrea.com/study/OpenFOAM/bc_settings.html)にて詳細に解説が為されているところではありますが、個人的にまとめておきたい情報をここに残しておこうと思います。

なお、本記事の環境は Foundation版(ver. 8) になります。ESI版とは若干入力方法が異なるところがありますのでご注意ください。

# 乱流パラメーターの算出
#### 乱流エネルギー
乱流エネルギー (turbulent kinetic energy)は以下の式で求められます。

```math
k = \frac{3}{2} (U I) ^ 2 
```

ここで_I_は乱流強度(turbulent intensity)であり、一般的なケース（十分に発達した流れ）であれば数パーセントです。ちなみにCFXなどではデフォルト値として 5% が設定されています。

#### 乱流エネルギー散逸率
乱流エネルギー消散速度(the rate of dissipation of turbulent kinetic energy)は次の式で算出できます。

```math 
\epsilon = \frac{C_{\mu}^{3/4} k^{3/2}}{l_m}
```

$C_{\mu}$ はモデルパラメーターであり、通常は **0.09**が使用されます。
$l_m$は混合長(mixing lenght)であり、管内流れで十分に流れが発達している場合は、

```math 
l_m = 0.07 L
```

で求められます。ここで *L* は代表長さであり、円管を対象とするときは直径を、矩形などを対象とする際は、$4 \times$ (面積) / (濡れぶち長さ(Wetted perimeter)) で算出します。

#### 乱流エネルギー比散逸率
$\omega$ は(the specific rate of dissipation of the turbulence kinetic energy k into internal thermal energy)は以下の式で求められます。

```math 
\omega = \frac{\epsilon}{C_{\mu} k} = \frac{k^{1/2}}{C_{\mu} ^{1/4} l_m}
```

# 乱流境界条件の設定例
以下ではOpenFOAMの入力の一例を示します。

#### _k_ 
```c++:k
inlet
{
   type     turbulentIntensityKineticEnergyInlet;
   intensity   0.05;  // 乱流強度を個々に入れる
   value   uniform 1; // $internalField;
}
outlet
{
   type    zeroGradient;
}
wall
{
   type    kqRWallFunction; // high Re flow
   // -
   type    kLowReWallFunction;// low - high Re flow
}
```
- [kqRWallFunction](https://cpp.openfoam.org/v8/classFoam_1_1kLowReWallFunctionFvPatchScalarField.html#details) は、高レイノルズ数流れの場合の壁関数モデル（k,q,R）になります（ゼロ勾配条件の単純なラッパーのようです）。
- [kLowReWallFunction](https://cpp.openfoam.org/v8/classFoam_1_1kLowReWallFunctionFvPatchScalarField.html#details)は、低レイノルズ数および高レイノルズ数乱流の場合の乱流運動エネルギー壁関数モデルです。
このモデルでは2つのモード、層流、乱流を、対応する nutWallFunction で指定された kappa と E から計算された y+ 値に基づいて切り替えています。

 
#### _[epsilon](https://cpp.openfoam.org/v8/classFoam_1_1turbulentMixingLengthDissipationRateInletFvPatchScalarField.html#details)_

```c++:epsilon
inlet
{
   type    turbulentMixingLengthDissipationRateInlet;
   mixingLength   0.001; // 混合長を入力する
   value   uniform 200;// $internalField;
}
outlet
{
   type    zeroGradient;
}
wall
{
   type    epsilonWallFunction;
}
```

- [epsilonWallFunction](https://cpp.openfoam.org/v8/classFoam_1_1epsilonWallFunctionFvPatchScalarField.html#details) は高Reおよび低Re乱流の壁面散逸モデルになります。
このモデルは、対応する `nutWallFunction` で指定された `kappa` と `E` から導出される `y+` 値に基づいて、層流と乱流を切り替えます。

#### _nut_

```
inlet
{
   type    calculated;
   value   uniform 0;
}
outlet
{
   type    zeroGradient;
   //-
   type    calculated;
   value   uniform 0;
}
wall
{
   type    nutWallFunction;
   value   uniform 0;
   // -
   type    nutkWallFunction;
   value   uniform 0;
}
```

- [nutWallFunction](https://cpp.openfoam.org/v8/classFoam_1_1nutWallFunctionFvPatchScalarField.html#details) は壁関数モデルになります。デフォルトyパラメーターは以下の通りです。

| Property | Description | Required | Default value |
| --- | --- | --- | --- |
| Cmu | $C_{mu}$ coefficient | no | 0.09 |
| kappa	| Von Karman constant | no | 0.41 |
| E | E coefficient | no | 9.8 |

これらパラメーターは境界条件設定で省略可能です（上記の例のようにnutWallFucntion を叩くのみでよいです）。

- [nutkWallFunction](https://cpp.openfoam.org/v8/classFoam_1_1nutkWallFunctionFvPatchScalarField.html#details) は高Re流れにおける壁関数モデルになります。

# さいごに
理論の部分などは省略してしまっています。順次、追記していく予定でいます。
わざわざプログラムにしなくてもという話もありますが、パラメーター計算用のpythonプログラムを[こちら](https://gist.github.com/greenty5/75b4ff0763bbf03aae48e26863d03601)にあげておきます。