OpenCFD社の OpenFOAM v2512 で 二層 $k-\varepsilon$ モデルが追加された。今後、役に立つかも知れないので、調査および動作確認を行う。
Turbulence: New two-layer wall treatment for kEpsilon model
1. 標準 k-ε モデル
最初に元となる、標準 $k-\varepsilon$ モデルの OpenFOAM での実装について記す。
1.1 ソースコード
ソースコードは次の場所に置かれている。
$FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kEpsilon
1.2 乱流動粘性係数
\nu_t = C_\mu \frac{k^2}{\varepsilon} \tag{1}
ここで、
$\nu_t :$ 乱流動粘性係数 $[m^2 s^{-1}]$
$C_{\mu} :$ 乱流動粘性係数のモデル係数 $[-]$
1.3 乱流運動エネルギー散逸率
\frac{\partial}{\partial t} \bigl( \alpha \rho \varepsilon \bigr) + \nabla \cdot \bigl( \alpha \rho \boldsymbol{u} \varepsilon \bigr) - \nabla^2 \bigl( \alpha \rho D_{\varepsilon} \varepsilon \bigr) = C_1 \alpha \rho \frac{G}{\nu_t} C_{\mu} k - \biggl( \frac{2}{3} C_1 - C_{3,RDT} \biggr) \alpha \rho \nabla \cdot \boldsymbol{u} \varepsilon - C_2 \alpha \rho \frac{\varepsilon}{k} \varepsilon + S_\varepsilon + S_{fvOptions} \tag{2}
D_{\varepsilon}=\frac{\nu_t}{\sigma_{\varepsilon}}+\nu \tag{3}
ここで、
$\alpha :$ 与えられた相の割合 $[-]$
$\rho :$ 流体の密度 $[kgm^{-3}]$
$G :$ レイノルズ応力テンソルの異方的成分による乱流運動エネルギー生成率 $[m^2 s^{-3}]$
$D_{\varepsilon} :$ $\varepsilon$ の有効拡散係数 $[m^2 s^{-1}]$
$\sigma_\varepsilon :$ モデル係数 $[-]$
$C_1 :$ モデル係数 $[s]$
$C_2 :$ モデル係数 $[-]$
$C_{3,RDT} :$ Rapid-distotion理論の圧縮項係数 $[-]$
$S_{\varepsilon} :$ $\varepsilon$ の内部ソースターム
$S_{fvOptions} :$ $fvOptions$ ディクショナリによるソースターム
1.4 乱流運動エネルギー方程式
\frac{\partial}{\partial t} \bigl( \alpha \rho k \bigr) + \nabla \cdot \bigl( \alpha \rho \boldsymbol{u} k \bigr) - \nabla^2 \bigl( \alpha \rho D_{k} k \bigr) = \alpha \rho G - \frac{2}{3} \alpha \rho \nabla \cdot \boldsymbol{u} k - \alpha \rho \frac{\varepsilon}{k} k + S_k + S_{fvOptions} \tag{4}
D_{k}=\frac{\nu_t}{\sigma_{k}}+\nu \tag{5}
ここで、
$D_k :$ $k$ の有効拡散係数 $[m^2 s^{-1}]$
$\sigma_k :$ モデル係数 $[-]$
$S_{k} :$ $k$ の内部ソースターム
$S_{fvOptions} :$ $fvOptions$ ディクショナリによるソースターム
注記
- 浮力の寄与は含まれていない
1.5 モデル係数のデフォルト値
C_\mu=0.09; C_1=1.44; C_2=1.92; C_{3,RDT}=0.0; \sigma_k=1.0; \sigma_{\varepsilon}=1.3
2. 二層 k-ε モデル
数値流体力学【第２版】 (p.91-92) に記載されている、二層 $k-\varepsilon$ モデルの説明を引用する。
二層モデルでは，低レイノルズ数での乱流流れにおける壁近傍の取扱いを改善し，記述する．これは，これまで説明した低レイノルズ数型 $k-\varepsilon$ モデルのように，粘性底層 ($y^+ < 1$) での壁近傍の格子点を配置し，壁に対して積分する．このモデルでは，壁に対する $k$ 方程式と，$\varepsilon$ 方程式の両方を積分するための低レイノルズ数型 $k-\varepsilon$ モデルに必要な非線形の壁面ダンピング関数に関連する数値安定性の問題 (Chen と Patel, 1998) を，境界層を次の二つの領域に分割することで回避した (Rodi, 1991)．
- 完全な乱流領域．$Re_y = y \sqrt k \nu \geq 200$：標準型 $k-\varepsilon$ モデルを用い、渦粘性を通常の関係式 $\mu_{t,t} = C_\mu \rho k^2 / \varepsilon$ で計算する。
- 粘性領域．$Re_y < 200$：この領域では，$k$ 方程式のみを解き，消散率 $\varepsilon = C_\mu^{3/4} k^{3/2} / l$ を評価するため，また，渦粘性 $\mu_{t,v} = C_\mu^{1/4} \rho \sqrt k l$，$A = 70$ を評価するために，$l = \kappa y [1-\exp(-Re_y/A)]$ を用いて，長さスケールを特定する．
(中略)
二層モデルは，初期の低レイノルズ数型 $k-\varepsilon$ モデルと比較して，格子依存性が減少し，数値安定性が向上している．また，流れの方程式の壁に対して積分する必要な (原文ママ) 複雑な流れの解析において，かなり人気の高いものとなっている．
(引用おわり)
二層 $k-\varepsilon$ モデルの OpenFOAM での実装について、以下に記す。
2.1 ソースコード
二層 $k-\varepsilon$ モデルのコードは、既存の $k-\varepsilon$ モデルに追加する形で実装されている。
2.2 モデルの説明
標準的な $k-\varepsilon$ モデルは高レイノルズ数モデルである。高レイノルズ数モデルであることは、たとえば $y^+ < 11$ である解像度が十分な計算グリッド上で、壁面せん断応力および熱流束の計算が正しくないことを意味する。二層 $k-\varepsilon$ モデルは、この問題に取り組むものである。壁に近い方を内層とし、乱流動粘性係数 $\nu_t$ および散逸率 $\varepsilon$ の代数的な関係式を定義する。一方で内層の外側と外層と定義し、$\nu_t$ および $\varepsilon$ は標準的な $k-\varepsilon$ モデルで計算する。
Two-layer wall treatment to allow operation on low-Re grids:
Jongen, T. and Marx, Y.P. (1997) “Design of an Unconditionally
Stable, Positive Scheme for the k−epsilon and Two-Layer
Turbulence Models”, Computers & Fluids, 26(5), 469-487.
2.3 モデルの係数
kEpsilon.H にモデルの係数のデフォルト値が記載されている。
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 0;
sigmak 1;
sigmaEps 1.3;
twoLayerTreatment false;
ReyStar 200.0;
ReyFactor 0.05;
}
下の3つが二層 $k-\varepsilon$ モデル向けに追加された変数である。
| 変数名 | 説明 |
|---|---|
| twoLayerTreatment | 二層 $k-\varepsilon$ モデル使用フラグ。 true ならば有効、false ならば無効である。デフォルト値は無効である。 |
| ReyStar | 内層と外層を分離するときに使用する乱流レイノルズ数。50.0 (低速流れ) から 200.0 (デフォルト値) の間で調整可能である。 |
| ReyFactor | 内層から外層への滑らかな移行。典型的には 0.05 (鋭い) から 0.2 (滑らか)。不安定なケースでは大きい値とすることが好ましい。 |
2.4 乱流動粘性係数
最初に式(1) の計算を行う。二層 $k-\varepsilon$ モデル使用フラグ (twoLayerTreatment) が無効 (false) であれば、計算はこれで終わりである。フラグが有効 (true) であれば $\nu_t$ の補正を行う。
現在のタイムステップの間での乱流レイノルズ数 $Re_y$ を固定するために乱流運動エネルギー (TKE) の前回値 $k_{old}$ を用いる。そして、その乱流レイノルズ数で内層 (粘性領域) と外層 (完全な乱流領域) を分離する。
Re_y = \biggr|\frac{y \sqrt{k_{old}}}{\nu} \biggl| \tag{6}
長さスケール $l_{\varepsilon}$ を次式で計算する。
l_\varepsilon = y Cl^{star} \biggr[1 - \exp \biggr( -\frac{Re_y}{2 Cl^{star}} \biggl) \biggl]\tag{7}
Cl^{star} = \kappa C_\mu^{-\frac{3}{4}} \tag{8}
ここで、$\kappa$ はカルマン定数で、値は $0.41$ である。
混合関数 $\lambda_{\varepsilon} = \lambda_{\varepsilon} (Re_y)$ を次式で計算する。
\lambda_\varepsilon = \frac{1}{2} \biggr[1 + \tanh \biggr( \frac{Re_y - Re_y^{star}}{A} \biggl) \biggl] \tag{9}
A = \frac{Re_y^{factor} Re_y^{star}}{\tanh^{-1} 0.98} \tag{10}
混合関数 $\lambda_{\varepsilon}$ は、次の特徴を持つ。
- 壁 $(y=0, Re_y=0)$ では $0$
- $Re_y \gg Re_y^{star}$ の場合、$1$
- $Re_y = Re_y^{star}$ 近傍で滑らかに遷移する
2.4.1 完全な乱流領域
式(1) および式(9) より、乱流動粘性係数 $\nu_t^{high}$ を計算する。
\nu_t^{high} = \lambda_\varepsilon \nu_t \tag{11}
2.4.2 粘性領域
粘性領域の長さスケールを計算して、乱流動粘性係数 $\nu_t^{low}$ を計算する。
A_\mu = 70\tag{12}
l_\mu = y Cl^{star} \biggr[1 - \exp \biggr(-\frac{Re_y}{A_\mu} \biggl) \biggl] \tag{13}
\nu_t^{low} = (1 - \lambda_\varepsilon) C_\mu l_\mu \sqrt{k} \tag{14}
式(11) および 式(14) の和を改めて乱流動粘性係数 $\nu_t$ とする。
\nu_t = \nu_t^{high} + \nu_t^{low} \tag{15}
2.5 乱流運動エネルギー方程式のソースターム
二層 $k-\varepsilon$ モデル使用フラグ (twoLayerTreatment) が有効 (true) であれば、次の量を $k$ 方程式のソースタームとして追加する。
S_k = - \alpha \rho (1 - \lambda_\varepsilon) \biggr(\frac{\sqrt k}{l_{\varepsilon}}-\frac{\varepsilon}{k} \biggl) k \tag{16}
式(4) に式(16) を代入すると、右辺の $-\alpha \rho \frac{\varepsilon}{k} k$ 項は次のようになる。
-\alpha \rho \frac{\varepsilon}{k} k \rightarrow -\alpha \rho \lambda_{\varepsilon} \frac{\varepsilon}{k} k -\alpha \rho (1 - \lambda_{\varepsilon}) \frac{\sqrt k}{l_{\varepsilon}} k \tag{17}
式(17) は、
- 完全な乱流領域の場合、$\lambda_\varepsilon=1$ で元々の乱流運動エネルギー方程式(式(4))に帰結する。
- 壁の場合、$\lambda_\varepsilon=0$ となり、$-\alpha \rho \frac{\varepsilon}{k} k$ 項は $-\alpha \rho \frac{\sqrt k}{l_{\varepsilon}} k$ となる。
3. 低レイノルズ数補正 (lowReCorrection)
最後に、$\varepsilon$ の壁関数のオプションである、低レイノルズ数補正について記す。
メッシュが粘性領域 $(y^+ < y_{lam}^+ \sim 11)$ で切られている場合、このオプションを有効にしないと $k-\varepsilon$ モデルで計算できないことがある。ここで、$y_{lam}^+$ は粘性領域と完全な乱流領域の境界である。lowReCorrection で有効・無効のコントロールを行う（デフォルトは無効）。
- lowReCorrection が有効の場合、段階的な混合が有効となる。
- lowReCorrection が無効の場合、完全な乱流領域のみ使用される。
3.1 ソースコード
ソースコードは次の場所に置かれている。
$FOAM_SRC/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions
3.1.1 y+ の計算
y^+ = C_{\mu}^{1/4} y \frac{\sqrt k}{\nu} \tag{18}
3.1.2 完全な乱流領域の壁面への ε の寄与
\varepsilon_{log} = C_{\mu}^{3/4} \frac{k^{3/2}}{\kappa y} \tag{19}
3.1.3 粘性領域の壁面への ε の寄与
\varepsilon_{vis} = 2 k \frac{\nu}{y^2} \tag{20}
3.1.4 段階的な ε の混合
lowReCorrection オプションが有効かつ $y^+ < y_{lam}^+$ の場合
\varepsilon_0 = \sum_{face} w_{face} \varepsilon_{vis, face} \tag{21}
その他の場合
\varepsilon_0 = \sum_{face} w_{face} \varepsilon_{log, face} \tag{22}
ここで、$w$ はコーナーウェイトである。
$\varepsilon_0$ は壁面での $\varepsilon$ がゼロ勾配条件を満たすように使用される。
4. 動作確認
OpenCFD社のチュートリアルケースを用いる。次のページを参照。
OpenFOAM のバージョンは、v2512 とする。
標準的な $k-\varepsilon$ モデルでは、$y^+<11$ のとき壁面せん断応力 $\tau$ の計算が正しく計算できないが、二層 $k-\varepsilon$ モデルではその計算が改善されることを確認する。
壁面せん断応力 $\tau$ より計算される摩擦係数 $C_f$で、壁面せん断応力の計算が改善されているかどうかを確認する。熱流束については今回使用する simpleFoam では取り扱えないため見送る。
4.1 チュートリアルケースの主な変更点
チュートリアルケースのディレクトリ構成を以下に示す
turbulentFlatPlate
├── Allclean
├── Allrun -> 4.1.1
├── README.md
├── plot -> 4.1.1
└── setups.orig
├── common
│ ├── 0.orig
│ │ ├── U
│ │ └── p
│ ├── Allclean
│ ├── Allrun
│ ├── Allrun-parallel
│ ├── Allrun.pre
│ ├── constant
│ │ └── transportProperties
│ └── system
│ ├── blockMeshDict
│ ├── blockMeshDict.template
│ ├── controlDict
│ ├── decomposeParDict
│ ├── fvSchemes
│ └── fvSolution -> 4.1.3
├── kEpsilon
│ ├── 0.orig
│ │ ├── epsilon
│ │ ├── k
│ │ └── nut
│ └── constant
│ └── turbulenceProperties
├── kEpsilonTwoLayer -> 4.1.2
│ ├── 0.orig
│ │ ├── epsilon
│ │ ├── k
│ │ └── nut
│ └── constant
│ └── turbulenceProperties
└── kOmegaSST
├── 0.orig
│ ├── k
│ ├── nut
│ └── omega
└── constant
└── turbulenceProperties
4.1.1 シェルスクリプト(Allrun および plot) の変更
- Allrun の 10 行目の kOmegaSST を kEpsilonTwoLayer に変更する。
- plot の 10 行目の kOmegaSST を kEpsilonTwoLayer に変更する。
4.1.2 二層 k-ε モデルの追加
- 上述の setups.orig 配下の kEpsilon ディレクトリをコピーして、kEpsilonTwoLayer ディレクトリを作成する。
- 二層 $k-\varepsilon$ モデルを使用するため、kEpsilonTwoLayer/constant/turbulenceProperties で twoLayerTreatment を true とする。
simulationType RAS;
RAS
{
RASModel kEpsilon;
turbulence on;
printCoeffs on;
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 0;
sigmak 1;
sigmaEps 1.3;
twoLayerTreatment true;
ReyStar 200.0;
ReyFactor 0.05;
}
}
- 低レイノルズ数補正 (lowReCorrection) の設定について、kEpsilonTwoLayer/0.orig/epsilon の bottomWall で使用する壁関数 epsilonWallFunction の lowReCorrection を on から off に変更する。
bottomWall
{
type epsilonWallFunction;
lowReCorrection off;
value $internalField;
}
4.1.3 代数方程式のソルバー設定 (system/fvSolutions)
- tolerance をすべて 1e-6 とする。
- residualControl をすべて 1e-4 とする。
- relaxationFactors を、p は 0.3、それ以外は 0.7 とする。
4.2 計算結果
各ケースの壁面摩擦係数 $C_f$ 対 x-レイノルズ数 $Re_x$ の計算結果を比較する。
$C_f$ および $Re_x$ の計算式は次の通り。
C_f = \frac{\sqrt{\tau_x^2 + \tau_y^2 + \tau_z^2}}{\frac{1}{2} U_{ref}^2} \tag{23}
Re_x = \frac{U_{ref} (x - x_0)}{\nu} \tag{24}
ここで、$U_{ref}=69.4 : [m/s]$、$x_0=0 :[m]$、$\nu = 1.388 \times 10^{-5} :[m^2/s]$ である。
比較する壁面摩擦係数の値は Weighardt の実験式に拠る。
C_{f,Weighardt}(Re_x) = 0.288 [log_{10}(Re_x)]^{-2.45} \tag{25}
4.2.1 標準 k-ε モデル
標準 $k-\varepsilon$ モデルの計算結果を以下に図示する。なお、$y^+=30, 50, 100$ のケースは、低レイノルズ数補正 (lowReCorrection) を無効にして計算した。$y^+=30, 50, 100$ のケースでは、$C_f$ は実験式と良く合うことがわかる。
4.2.2 二層 k-ε モデル
二層 $k-\varepsilon$ モデルの計算結果を以下に図示する。$y^+ < 11$ のケースで、$y^+ = 0.05$ のケースのみ $C_f$ が実験式と合わない（メッシュの高さがイマイチなのだろうか）が、他は実験式とよく合っていることがわかる。
5. まとめ
4.2節で図示した通り、標準的な $k-\varepsilon$ モデルと比較して、二層 $k-\varepsilon$ モデルでは、$y^+<11$ のときの壁面摩擦係数 $C_f$ (即ち、壁面せん断応力 $\tau$) の計算が改善されることを確認できた。