大学院生のための数値流体力学③：LES

SGS応力のモデル化の基礎となっているエネルギーカスケードという考え方について説明する

エネルギーカスケード

流れ場に存在する乱れ（渦）の生成と消滅の過程は次のようになる

1. 流れ場のスケール（壁面の段差や航空機の翼）からエネルギーが与えられて大きな渦が生成される
2. 大きな渦が分裂を繰り返して小さな渦になる
3. 最終的にこれ以上分裂できない最小スケールになった渦が粘性によって熱となって散逸する

このような，外部から大きい渦へ注入されたエネルギーが小さな渦へと段階的に伝達していく過程のことを エネルギーカスケード（Energy cascade）という

最初に作られた大きな渦は流れ場の影響を受けており，統計的な非等方性が現れる（大きな渦の形状や数は流れ場に依存する）が，渦が分裂するにしたがって最初の渦が持っていた非等方性はなくなっていき，何度も分裂を繰り返したある程度の小さいスケールの渦は流れ場に依存しない等方的な渦になると考えられる

このような考え方を局所等方性（局所＝小さいスケールでの）といい，流れ場の影響を受けない小さいスケールの渦はどんな流れ場に対しても普遍的なものであるという仮定を「コルモゴロフの第一相似則」という

局所等方性をもつような小さいスケールの渦は，流れ場によって作られる大きな渦と比べて時間スケール（分裂や合体にかかる時間）が非常に小さく，常に平衡状態にあるとみなせる

よって，局所等方性が成り立ち平衡状態にある普遍的な渦のスケールを「普遍平衡領域」とよぶ

また，流れ場に存在する乱れのエネルギーに注目すると，普遍平衡領域の渦のもつエネルギーはほとんどなく，その大部分は流れ場の形状に依存する大きな渦が担っている

普遍平衡領域に対して，流れ場の形状に依存し，乱れエネルギーの大部分を保有する大きな渦のスケールを「エネルギー保有領域」とよぶ

LESは，普遍平衡領域の乱れをモデル化し，エネルギー保有領域の乱れを直接格子上で計算することで，どんな流れ場にも使える普遍的なモデルとなっている

慣性小領域

渦が粘性によって熱として散逸していく渦のスケールを「散逸領域」といい，普遍平衡領域の中でも最も小さいスケールである

前々回の記事で説明したように，流れ場の代表スケール（エネルギー保有領域のスケール）\(l\)と流れ場に存在する渦の最小スケール（散逸領域のスケール）\(\eta\)には次のような関係がある

\begin{equation}
\frac{l}{\eta}={Re}^{\frac{3}{4}}
\end{equation}

この式は，レイノルズ数が大きくなるとその3/4乗に比例してエネルギー保有領域と散逸領域のスケールが離れていくことを意味している

レイノルズ数が十分に高い流れ場では，エネルギー保有領域と散逸領域には大きな隔たりがあり，大きな渦が分裂を繰り返して普遍平衡領域に入ってからしばらくはひたすら分裂を続けるだけの領域がある

このように，普遍平衡領域の中で粘性による散逸の影響を受けずに分裂だけを繰り返す領域を「慣性小領域」とよぶ

また，このような領域が存在するという仮説を「コルモゴロフの第二相似則」という

エネルギースペクトル

これまでの説明をグラフにしたものがエネルギースペクトルである

グラフの縦軸に乱れの速度をフーリエ変換して得られるエネルギースペクトル\(E(k)\)をとり，横軸に乱れの波数\(k\)をとると，次のようになる

グラフ左側の波数が小さい（波長が大きい）領域がエネルギー保有領域であり，流れ場からエネルギーをもらうとともに流れ場の大部分のエネルギーを保有している

グラフ右側の波数が大きい（波長が小さい）領域が散逸領域であり，粘性による散逸によってエネルギースペクトルが急速に減少している

グラフ中央の領域が慣性小領域であり，慣性小領域におけるエネルギースペクトルは次元解析によって次のように求められる

\begin{equation}
E(k)=\alpha \varepsilon^{2}{3}k^{-\frac{5}{3}}
\end{equation}

慣性小領域における-5/3乗スペクトルは実験でも確認されており，係数\(\alpha\)は実験的に \(\alpha=1.2 \sim 2.0\)であることがわかっている

LESのモデルにおいては，格子サイズを慣性小領域に設定することを前提としており，慣性小領域を境にGS成分とSGS成分が隔てられている

順カスケードと逆カスケード

乱れのエネルギーが渦の分裂によってより小さい渦へと伝達していく過程を「順カスケード」といい，その逆に，小さい渦が合体することによって乱れのエネルギーが大きな渦へと伝達していく過程を「逆カスケード」という

乱れの中では常に順カスケードと逆カスケードが起こっているが，それらは相対的に非常に小さな時間スケールで起こっており平衡状態にある

全体としては順カスケードのほうが多いため，平均的に見ると，渦はより小さい渦へと分裂していく順カスケードにあるとみなせる

SGS応力の式に \(u_{i}=\overline{u}_{i}+u'_{i}\) を代入すると，SGS応力は異なる役割を持つ3つの項に分けられる

\begin{align}
\tau_{ij}
&=\overline{\left(\overline{u}_{i}+u'_{i}\right)\left(\overline{u}_{i}+u'_{i}\right)}-\overline{u}_{i}\overline{u}_{j} \\
&=\overline{\overline{u}_{i}\overline{u}_{j}+\overline{u}_{i}u'_{j}+u'_{i}\overline{u}_{j}+u'_{i}u'_{j}}-\overline{u}_{i}\overline{u}_{j} \\
&=\left(\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{u}_{i}\overline{u}_{j}\right)+\left(\overline{\overline{u}_{i}u'_{j}}+\overline{u'_{i}\overline{u}_{j}}\right)+\overline{u'_{i}u'_{j}} \\
&=L_{ij}+C_{ij}+R_{ij}
\end{align}

（レナード項：Leonard Term）

\begin{equation}
L_{ij}=\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{u}_{i}\overline{u}_{j}
\end{equation}

レナード項は，平均的に見るとGS成分からSGS成分への運動エネルギーの輸送（順カスケード）を担っており，既知の変数\(\overline{u}_{i}\)にフィルターをかけることで直接計算することができる

（クロス項：crossTerm）

\begin{equation}
C_{ij}=\overline{\overline{u}_{i}u'_{j}}+\overline{u'_{i}\overline{u}_{j}}
\end{equation}

クロス項は，平均的に見るとSGS成分からGS成分への運動エネルギーの逆輸送（逆カスケード）を担っている

（SGSレイノルズ応力：SGS Reynolds stress）

\begin{equation}
R_{ij}=\overline{u'_{i}u'_{j}}
\end{equation}

SGSレイノルズ応力は，主にSGSにおける運動エネルギーの散逸を担っている

流れ場にフィルターをかけると，流れ全体の運動エネルギー \(\overline{k}=\frac{\overline{u_{i}u_{i}}}{2}\)（乱流エネルギーではない）はGS成分 \(k_{_{GS}}=\frac{\overline{u}_{i}\overline{u}_{i}}{2}\)とSGS成分 \(k_{_{SGS}}=\frac{\overline{u_{i}u_{i}}-\overline{u}_{i}\overline{u}_{i}}{2}\) に分けられる

運動エネルギーのそれぞれの成分の輸送方程式は次のように表される

\begin{align}
&\frac{\partial k_{_{GS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{GS}}}{\partial x_{j}}
=
\tau_{ij}\overline{S}_{ij}-\varepsilon_{_{GS}}
+\frac{\partial}{\partial x_{j}}
\left(-\overline{u}_{i}\tau_{ij}-\frac{\overline{p}\overline{u}_{j}}{\rho}+\nu\frac{\partial k_{_{GS}}}{\partial x_{j}}\right) \\
\\
&\frac{\partial k_{_{SGS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{SGS}}}{\partial x_{j}}
=
-\tau_{ij}\overline{S}_{ij}-\varepsilon_{_{SGS}} \\
&\qquad +\frac{\partial}{\partial x_{j}}
\left\{\overline{u}_{i}\tau_{ij}-\frac{\overline{u_{i}u_{i}u_{j}}+\overline{u}_{j}\overline{u_{i}u_{i}}}{2}-\frac{\overline{pu_{j}}-\overline{p}\overline{u}_{j}}{\rho}+\nu\frac{\partial k_{_{SGS}}}{\partial x_{j}}\right\}
\end{align}

GS成分の輸送方程式においては\(\tau_{ij}S_{ij}\)の項がSGS運動エネルギーへの輸送率を表し，これを受けてSGS成分の輸送方程式では\(-\tau_{ij}S_{ij}\)の項がSGS運動エネルギーの生成率になる

SGSにおいて運動エネルギーの生成と散逸がつりあっているという局所平衡を仮定すると，SGSにおける運動エネルギーの散逸率\(\varepsilon_{_{SGS}}\)は次のようになる

\begin{equation}
\varepsilon_{_{SGS}}=-\tau_{ij} \overline{S}_{ij}
\tag{i}
\end{equation}

ここで突然SGS応力の渦粘性近似を行う

\begin{equation}
\tau_{ij}=\frac{2}{3}k\delta_{ij}-2\nu_{_{SGS}}\overline{S}_{ij}
\tag{ii}
\end{equation}

(i)に(ii)を代入すると次のようになる

\begin{align}
\varepsilon_{_{SGS}}
=-\left(\frac{2}{3}k\delta_{ij}-2\nu_{_{SGS}}\overline{S}_{ij}\right) \overline{S}_{ij}
=2\nu_{_{SGS}}\overline{S}_{ij}\overline{S}_{ij}-\frac{2}{3}k\delta_{ij}\overline{S}_{ij}
=2\nu_{_{SGS}}\overline{S}_{ij}\overline{S}_{ij}
\end{align}

（\(\delta_{ij}\overline{S}_{ij}=\overline{S}_{kk}=0\)）

SGS渦粘性係数\(\nu_{_{SGS}}\)は速度\(q\)と長さ\(l\)の次元 \(\mathrm{[m^2/s]}\) を持つので，次元解析的に次のように表すことができる

\begin{equation}
\nu_{_{SGS}}=C_{\nu}ql
\tag{iii}
\end{equation}

同様に，エネルギー散逸率の次元 \(\mathrm{[m^2/s^3]}\) を考慮すると，\(\varepsilon_{_{SGS}}\)は次元解析的に次のように表せる

\begin{equation}
\varepsilon_{_{SGS}}=2\nu_{_{SGS}}\overline{S}_{ij}\overline{S}_{ij}=C_{\varepsilon}\frac{q^{3}}{l}
\tag{iv}
\end{equation}

(iii)と(iv)より

\begin{equation}
\left\{\begin{array}{rcl}
q&=&\nu_{_{SGS}}/C_{\nu}l \\
q&=&\left(\frac{1}{C_{\varepsilon}}2\nu_{_{SGS}}l \overline{S}_{ij}\overline{S}_{ij}\right)^{\frac{1}{3}}
\end{array}\right.
\end{equation}

\begin{align}
\frac{\nu_{_{SGS}}}{C_{\nu}l}=\left(\frac{1}{C_{\varepsilon}}2\nu_{_{SGS}}l \overline{S}_{ij}\overline{S}_{ij}\right)^{\frac{1}{3}}
\rightarrow \quad
\nu_{_{SGS}}=\left\{\left(\frac{C_{\nu}^{3}}{C_{\varepsilon}}\right)^{\frac{1}{4}}l\right\}^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}
\end{align}

上式に長さのスケール\(l\)としてフィルター幅\(\Delta\)を用いて，定数をまとめると，Smagorinskyモデルの\(\nu_{_{SGS}}\)を求めることができる

\begin{equation}
\nu_{_{SGS}}=\left(C_{s}\Delta\right)^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}
\tag{v}
\end{equation}

Smagorinsky定数の導出

(v)を(iv)に代入すると，次のようになる

\begin{align}
&\varepsilon_{_{SGS}}=2\left\{\left(C_{s}\Delta\right)^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}\right\}\overline{S}_{ij}\overline{S}_{ij}
=\left(C_{s}\Delta\right)^{2}\left(2\overline{S}_{ij}\overline{S}_{ij}\right)^{\frac{3}{2}} \\
\\ &\rightarrow\quad
|\mathbf{S}|^{2}=\overline{S}_{ij}\overline{S}_{ij}=\frac{\varepsilon_{_{SGS}}^{2/3}}{2\left(C_{s}\Delta\right)^{4/3}}
\tag{vi}
\end{align}

ここで，フィルター幅\(\Delta\)が慣性小領域にあるように格子分割を定めると，その波数\(k\)付近のエネルギースペクトル\(E(k)\)は-5/3乗則 \(E(k)=\alpha\varepsilon^{\frac{2}{3}}k^{-\frac{5}{3}}\) に従うと考えられる（\(k\)は波数で\(\varepsilon\)は粘性散逸総量）

このとき，\(|\mathbf{S}|^{2}=\overline{S}_{ij}\overline{S}_{ij}\) は渦度スペクトル関数の波数積分から次のように求められる

\begin{align}
|\mathbf{S}|^{2}=\int_{0}^{\pi/4} k^{2}E(k)\:dk
=\frac{3}{4}\alpha|\varepsilon|^{\frac{2}{3}}\left(\frac{\pi}{\Delta}\right)^{\frac{4}{3}}
\tag{vii}
\end{align}

SGSスケールの運動エネルギー散逸率とエネルギースペクトルにおける粘性散逸総量が等しい \(\varepsilon_{_{SGS}}=\varepsilon\) とすると，(vi)と(vii)よりSmagorinsky定数が計算できる

\begin{align}
\frac{\varepsilon_{_{SGS}}^{2/3}}{2\left(C_{s}\Delta\right)^{4/3}}
=\frac{3}{4}\alpha|\varepsilon_{_{SGS}}|^{\frac{2}{3}}\left(\frac{\pi}{\Delta}\right)^{\frac{4}{3}}
\quad \rightarrow \quad
C_{s}=\frac{1}{\pi}\left(\frac{2}{3\alpha}\right)^{\frac{3}{4}}
\end{align}

ここで，実験的に \(\alpha=1.5\) であることが確認されているので，Smagorinsky定数の理論値 \(C_{s}=0.173\) が導かれる

参考
≫Smagorinsky, J.: General circulation experiments with the primitive equations. I. The basic experiment. Mon. Weather Rev. 91, 99–164 (1963)

エネルギーカスケードの考え方より，SGSの流れの挙動は，GSの流れの中でもより小さく，SGSに近いスケールの流れによって決定されていると考えられる

そこで，GSの中からより小さいスケールを取り出すため，Filtered Naver-Stokes方程式にさらにフィルターをかけることを考える

\begin{align}
\frac{\partial \widetilde{\overline{u}}_{i}}{\partial t}
+\frac{\partial \widetilde{\overline{u_{i}u_{j}}}}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \widetilde{\overline{p}}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{2\nu\widetilde{\overline{S}}_{ij}\right\}
\end{align}

両辺から\(\frac{\partial \widetilde{\overline{u_{i}u_{j}}}}{\partial x_{j}}\)を引いて\(\frac{\partial \widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}}{\partial x_{j}}\)を足すと

\begin{align}
\frac{\partial \widetilde{\overline{u}}_{i}}{\partial t}
+\frac{\partial \widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \widetilde{\overline{p}}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{2\nu\widetilde{\overline{S}}_{ij}-T_{ij}\right\}
\end{align}

サブテストフィルタースケール応力\(T_{ij}\)（GSの中の小さいスケールの乱れ）

\begin{equation}
T_{ij}=\widetilde{\overline{u_{i}u_{j}}}-\widetilde{\overline{u}_{i}\overline{u}_{j}}
\end{equation}

が求まる

サブテストフィルタースケール応力がSGSの乱れに与える影響を考えるために，\(T_{ij}\)と\(\tau_{ij}\)の関係式を求める

\(\tau_{ij}\)にテストフィルターをかけると

\begin{align}
\widetilde{\tau}_{ij}=\widetilde{\overline{u_{i}u_{j}}}-\widetilde{\overline{u}_{i}\overline{u}_{j}}
=\widetilde{\overline{u_{i}u_{j}}}-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}
-\left(\widetilde{\overline{u}_{i}\overline{u}_{j}}-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}\right)
=T_{ij}-L_{ij}
\end{align}

\(T_{ij}\)と\(\tau_{ij}\)の関係式

\begin{equation}
L_{ij}=\widetilde{\overline{u}_{i}\overline{u}_{j}}-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}=T_{ij}-\widetilde{\tau}_{ij}
\tag{i}
\end{equation}

を求めることができる

\(L_{ij}\)は，サブテストフィルタースケール応力に対するレナード項にあたり，既知の速度成分\(\overline{u}_{i}\)のみで構成されているので直接計算で求めることができる（\(\overline{u_{i}u_{j}}\)は二次の相関項なので計算できない）

\(T_{ij}\)と\(\tau_{ij}\)の関係式からSmagorinsky定数を計算するために，それぞれのスケールの乱流応力にSmagorinskyモデルを適用する

\begin{align}
T_{ij}-\frac{1}{3}T_{kk}\delta_{ij}
&=-2C\tilde{\Delta}^{2}|\widetilde{\overline{S}}|\widetilde{\overline{S}}_{ij} \\
\tilde{\tau}_{ij}-\frac{1}{3}\tilde{\tau}_{kk}\delta_{ij}
&=-2C\overline{\Delta}^{2}\widetilde{|\overline{S}|\overline{S}_{ij}}
\tag{ii}
\end{align}

\(\tilde{\Delta}\)はテストフィルターのフィルター幅，\(\overline{\Delta}\)はグリッドフィルターのフィルター幅である

(ii)を(i)に代入すると次のようになる

\begin{equation}
L_{ij}-\frac{1}{3}\left(T_{kk}-\tilde{\tau}_{kk}\right)\delta_{ij}
=T_{ij}-\frac{1}{3}T_{kk}\delta_{ij}-\left(\widetilde{\tau}_{ij}-\frac{1}{3}\widetilde{\tau}_{kk}\delta_{ij}\right) \\
\end{equation}

\begin{equation}
\rightarrow \quad
L_{ij}^{a}
=-2C\tilde{\Delta}^{2}|\widetilde{\overline{S}}|\widetilde{\overline{S}}_{ij}-\left(-2C\overline{\Delta}^{2}\widetilde{|\overline{S}|\overline{S}_{ij}}\right)
=-2C\overline{\Delta}^{2}\left(\frac{\tilde{\Delta}^{2}}{\overline{\Delta}^{2}}|\widetilde{\overline{S}}|\widetilde{\overline{S}}_{ij}-\widetilde{|\overline{S}|\overline{S}_{ij}}\right)
\end{equation}

\(L_{ij}^{a}=L_{ij}-\frac{1}{3}L_{kk}\delta_{ij}\) は \(L_{ij}\) の非等方テンソル（anisotropic tensor）であり，テストフィルターとグリッドフィルターのフィルター幅の比を \(\gamma=\tilde{\Delta}/{\overline{\Delta}}\)とすれば次のようになる

\begin{equation}
L_{ij}^{a}=-2C\overline{\Delta}^{2}M_{ij}
\tag{iii}
\end{equation}

\begin{equation}
M_{ij}=\gamma^{2}|\widetilde{\overline{S}}|\widetilde{\overline{S}}_{ij}-\widetilde{|\overline{S}|\overline{S}_{ij}}
\end{equation}

ここで，\(L_{ij}^{a}\) と \(-2C\overline{\Delta}^{2}M_{ij}\) は対称テンソルかつトレースレス \(L_{kk}=-2C\overline{\Delta}^{2}M_{kk}=0\) なので，それぞれの独立な成分は5つである

１つのスカラー\(C\)だけを用いて２つのテンソルの独立な５つの成分をすべて一致させることはできないので，２つのテンソルの誤差が最小になるように，最小二乗法を用いて\(C\)を決定する

(iii)の両辺を比較したときの誤差テンソルを \(e_{ij}=L_{ij}^{a}-\left(-2\overline{\Delta}^{2}M_{ij}\right)\) としたとき，最小二乗法より\(e_{ij}\)の二乗を考えると

\begin{align}
&E=e_{ij}e_{ij}=\left(L_{ij}^{a}+2C\overline{\Delta}^{2}M_{ij}\right)^{2}\\
\\
\rightarrow \quad
&\frac{\partial E}{\partial C}
=2\left(L_{ij}^{a}+2C\overline{\Delta}^{2}M_{ij}\right)2\overline{\Delta}^{2}M_{ij}
=4\overline{\Delta}^{2}\left(L_{ij}^{a}M_{ij}+2C\overline{\Delta}^{2}M_{ij}M_{ij}\right)
\end{align}

となり，誤差を最小にするには \(\frac{\partial E}{\partial C}=0\) となればいいので

\begin{align}
&\frac{\partial E}{\partial C}
=4\overline{\Delta}^{2}\left(L_{ij}^{a}M_{ij}+2C\overline{\Delta}^{2}M_{ij}M_{ij}\right)=0 \\
\\ \rightarrow \quad
&L_{ij}^{a}M_{ij}+2C\overline{\Delta}^{2}M_{ij}M_{ij}=0
\end{align}

よって，Dynamic Smagorinskyモデルの係数\(C\)は次のようになる

\begin{equation}
C=-\frac{1}{2\overline{\Delta}^{2}}\frac{L_{ij}^{a}M_{ij}}{M_{ij}M_{ij}}
\end{equation}

実際の計算では\(C\)の値の変動を緩和するため，右辺の分母分子にそれぞれ何らかの平均値を用いる必要がある

\begin{equation}
C=-\frac{1}{2\overline{\Delta}^{2}}\frac{\left[L_{ij}^{a}M_{ij}\right]}{\left[M_{ij}M_{ij}\right]}
\end{equation}

参考
≫Lilly, D.K. A proposed modification of the Germano subgrid-scale closure method (1992) Physics of Fluids A, 4 (3), pp. 633-635.

Dynamic Smagorinskyモデルと同様に，サブテストフィルタースケール応力のレナード項\(L_{ij}\)は次のようになる

\begin{equation}
L_{ij}
=\widetilde{\overline{u}_{i}\overline{u}_{j}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}
\end{equation}

ここで，テストフィルターをかけた流れ場における流体の運動エネルギー\(\widetilde{k}\)とエネルギー散逸率\(\widetilde{\varepsilon}\)を考えると，次のようにテストフィルタースケールとサブテストフィルタースケールに分けることができる

\begin{align}
&\widetilde{k}=\frac{1}{2}\widetilde{\overline{u}_{i}\overline{u}_{i}}
=\frac{1}{2}\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{i}
+\frac{1}{2}\left(\widetilde{\overline{u}_{i}\overline{u}_{i}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{i}\right)
\\ \\
&\widetilde{\varepsilon}=\left(\nu+\nu_{_{SGS}}\right)\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}
=\left(\nu+\nu_{_{SGS}}\right)\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}+\left(\nu+\nu_{_{SGS}}\right)\left(\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}
-\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\right)
\end{align}

上式の中で，サブテストフィルタースケール（GSの中の小さいスケール）の運動エネルギー\(k_{_{test}}\)とエネルギー散逸率\(\varepsilon_{_{test}}\)は次の部分である

\begin{align}
&k_{_{test}}=\frac{1}{2}\left(\widetilde{\overline{u}_{i}\overline{u}_{i}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{i}\right)=\frac{L_{ii}}{2}
\\ \\
&\varepsilon_{_{test}}=
\left(\nu+\nu_{_{SGS}}\right)\left(\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}
-\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\right)
\tag{i}
\end{align}

サブテストフィルタースケールは壁面から離れた十分にレイノルズ数の高い領域なので，エネルギーの散逸は分子粘性と渦粘性の両方によって担われている（\(\nu+\nu_{_{SGS}}\)）

\(C_{k}\)

サブテストフィルタースケール応力のレナード項\(L_{ij}\)とSGS応力\(\tau_{ij}\)は相似の関係にあることがわかっている

そこで，\(L_{ij}\)に対してk Equationモデルを適用する

\begin{align}
&L_{ij}-\frac{1}{3}L_{kk}\delta_{ij}
=-2C_{k}\widetilde{\Delta}\sqrt{k_{_{test}}}\widetilde{\overline{S}}_{ij} \\
\rightarrow \quad
&L_{ij}^{a}=-2C_{k}M_{ij}, \qquad
M_{ij}=\widetilde{\Delta}\sqrt{k_{_{test}}}\widetilde{\overline{S}}_{ij}
\end{align}

この式に対して，Dynamic Smagorinskyモデルと同様に最小二乗法を用いると

\begin{equation}
e_{ij}=L_{ij}^{a}-\left(-2C_{k}M_{ij}\right), \quad
E=e_{ij}e_{ij}=\left(L_{ij}^{a}+2C_{k}M_{ij}\right)^{2}
\end{equation}

\begin{equation}
\frac{\partial E}{\partial C_{k}}
=2\left(L_{ij}^{a}+2C_{k}M_{ij}\right)2M_{ij}
=4\left(L_{ij}^{a}M_{ij}+2C_{k}M_{ij}M_{ij}\right)=0 \\
\end{equation}

\(C_{k}\)は次のように導かれる

\begin{equation}
C_{k}=-\frac{1}{2}\frac{L_{ij}M_{ij}}{M_{kl}M_{kl}}
\end{equation}

\(C_{\varepsilon}\)

\(\varepsilon_{_{SGS}}\)と同様に，\(\varepsilon_{_{test}}\)は次の式で表される

\begin{equation}
\varepsilon_{_{test}}=C_{\varepsilon}\frac{k_{_{test}}^{\frac{3}{2}}}{\widetilde{\Delta}}
\tag{ii}
\end{equation}

(i)と(ii)より，\(C_{\varepsilon}\)を導くことができる

\begin{align}
C_{\varepsilon}\frac{k_{_{test}}^{\frac{3}{2}}}{\widetilde{\Delta}}
=
\left(\nu+\nu_{_{SGS}}\right)\left(\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}
-\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\right) \\
\\ \rightarrow \quad
C_{\varepsilon}
=
\frac{\left(\nu+\nu_{_{SGS}}\right)\widetilde{\Delta}}{k_{_{test}}^{3/2}}\left(\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}
-\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\right) \\
\end{align}

参考
≫OpenFOAM User Guide Dynamic k eqn
≫W-W Kim and S. Menon. A new dynamic one-equation subgrid-scale model for large eddy simulations. In 33rd Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 1995.

以下のサイトにあるPDFのpp.184-193が上記の論文
≫A New Approach to Validate Subgrid Models in Complex High Reynolds Number Flows.

スケール相似則

SGS成分の定義式 \(u'_{i}=u_{i}-\overline{u}_{i}\) の両辺にさらにフィルターをかけると次のようになる

\begin{equation}
\overline{u'}_{i}=\overline{u}_{i}-\overline{\overline{u}}_{i}
\tag{i}
\end{equation}

左辺の\(\overline{u'}_{i}\)は，フィルターをかけることによって粗視化された変動量であり，\(u'_{i}\)のうちの大きなスケールを意味する

右辺第2項の\(\overline{\overline{u}}_{i}\)はフィルターをかけることによってさらに粗視化されたGS成分であり，\(\overline{u}_{i}\)のうちのより大きなスケールを意味しているので，右辺全体の \(\overline{u}_{i}-\overline{\overline{u}}_{i}\) はGS成分とGS成分のうちのより大きなスケールとの差，すなわちGS成分のうちの小さいスケールを意味している

以上の考察から，(i)の等式は，SGS成分\(u'_{i}\)のうちの大きなスケールとGS成分\(\overline{u}_{i}\)のうちの小さいスケールが相似であることを意味しており，このことをスケール相似則という

Bardinaモデル

上記のスケール相似則の考え方を用いて，BardinaはSGS応力のクロス項とSGSレイノルズ応力項にたいして以下のモデルを提案した

（クロス項：crossTerm）

\begin{equation}
C_{ij}=\overline{\overline{u}_{i}u'_{j}}+\overline{u'_{i}\overline{u}_{j}}
=\overline{\overline{u}}_{i}\left(\overline{u}_{j}-\overline{\overline{u}}_{j}\right)
+\left(\overline{u}_{i}-\overline{\overline{u}}_{i}\right)\overline{\overline{u}}_{j}
\end{equation}

（SGSレイノルズ応力：SGS Reynolds stress）

\begin{equation}
R_{ij}=\overline{u'_{i}u'_{j}}
=C_{B}\left(\overline{u}_{i}-\overline{\overline{u}}_{i}\right)\left(\overline{u}_{j}-\overline{\overline{u}}_{j}\right), \qquad C_{B}=1
\end{equation}

ここで，レナード項

\begin{equation}
L_{ij}=\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{u}_{i}\overline{u}_{j}
\end{equation}

はGS成分の速度にフィルターをかけることで直接計算できるので，Bardinaによるスケール相似則モデルは次のようになる

\begin{align}
\tau_{ij}
&=L_{ij}+C_{ij}+R_{ij} \\
&=\left(\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{u}_{i}\overline{u}_{j}\right)
+\left\{\overline{\overline{u}}_{i}\left(\overline{u}_{j}-\overline{\overline{u}}_{j}\right)+\left(\overline{u}_{i}-\overline{\overline{u}}_{i}\right)\overline{\overline{u}}_{j}\right\}
+\left\{\left(\overline{u}_{i}-\overline{\overline{u}}_{i}\right)\left(\overline{u}_{j}-\overline{\overline{u}}_{j}\right)\right\} \\
&=\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{\overline{u}}_{i}\overline{\overline{u}}_{j}
\end{align}