第13-2回　Large eddy simulationにおける圧縮性NS方程式[julia]

*第13-5回にコードがあります。

今回はLESにおける支配方程式を導出します。

添え字や用語として下記の二つを使用します。

GS : Grid scale
SGS : Sub-Grid scale

GSは計算格子サイズ以上を指し、SGSは計算格子サイズ未満を指します。

まず、圧縮性NS方程式は下記の通りです。質量保存則、運動量保存則、運動エネルギー保存則の順に記載しています。

\begin{equation}
\frac{\partial \rho}{\partial t} + \frac{\partial }{\partial x_j}\left(\rho u_j\right) = 0
\end{equation}

\begin{equation}
\frac{\partial \left(\rho u_j\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_i u_j + \delta_{ij}p \right) = \frac{\partial \sigma_{ij}}{\partial x_j}
\end{equation}

\begin{equation}
\frac{\partial e}{\partial t} + \frac{\partial }{\partial x_j}\left[\left(e+ p\right)u_j\right] =
\frac{\partial }{\partial x_j}\left( \sigma_{ij} u_{j}\right)
-\frac{\partial }{\partial x_j}q_{j}
\end{equation}

ただし、下記の通りの定義です。

\begin{equation}
e = \frac{p}{\gamma – 1} + \frac{1}{2} \rho u_{i} u_{i}
\end{equation}

\begin{equation}
p = \rho R T
\end{equation}

\begin{equation}
\sigma_{ij} = 2 \mu S_{ij} – \frac{2}{3} \mu \delta_{ij} S_{kk}
\end{equation}

\begin{equation}
S_{ij} = \frac{1}{2} \left( \frac{\partial u_{j}}{\partial x{i}} + \frac{\partial u_{i}}{\partial x{j}}\right)
\end{equation}

\begin{equation}
q_{j} = \lambda \frac{\partial T}{\partial x_j}
\end{equation}

LESの特徴はSGS以下はモデル化し、GSの渦は直接解くことにあります。

そのため、GSとSGSを分離する必要があります。

そこで下記のようにフィルタリング関数(G)との畳み込み積分を行います。

\begin{eqnarray}
\overline{u} = \int^{\infty} _{-\infty} G(x_i-x’) u(x’) dx’
\end{eqnarray}

省略のため一つの積分を書きましたが、本来は三方向の積分を行います。

フィルタリング関数にもよりますが、セル中心から一定以上離れた場合すなわち\(x_i-x’\)が0からずれるほど\(G=0\)に近づく構造になっています。

バーがついているものはGSを指します。

ただし，このフィルタリング関数はGSとSGSを区別する概念だと思ってください．

今回使用するモデルでは，具体的なフィルタリング関数は使いません．

質量保存則に対して、フィルタリング関数を適用してみます。積分と微分の順序は入れ替えても問題ありません。（証明は省きます。）

\begin{equation}
\overline{\frac{\partial \rho}{\partial t}} + \overline{\frac{\partial }{\partial x_j}\left(\rho u_j\right)} = 0
\nonumber
\end{equation}

\begin{equation}
\Leftrightarrow
\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial }{\partial x_j}\left(\overline{\rho u_j}\right) = 0
\end{equation}

左辺第一項はバーが付いているが、セルサイズ程度の値であるため、コードの中では今まで通りセル中心の値として使用します。つまりコード上では何も変化がありません。

一方、左辺第二項は積であり、このままでは評価することができません。

そのため、次のファーブル平均(Favre Average,favre-filtering)を導入します．

これは，質量荷重平均であり，次の通り定義されます．

\begin{eqnarray}
\tilde{\phi} = \frac{\overline{\rho \phi}}{\bar{\rho}}
\end{eqnarray}

これを導入する理由は，「余計な項を増やさないから」，「RANSの手法を応用できるから」です．

もし \(\overline{\rho u_j} = \bar{\rho} \bar{u_j} + \alpha\) とすると余計なSGS項が増えてしまいます．

以上のファーブル平均を質量保存則の左辺第2項に適用すると次のようになります．

\begin{equation}
\Leftrightarrow
\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial }{\partial x_j}\left(\overline{\rho} \tilde{u_j}\right) = 0
\end{equation}

チルダがついていますが、コードの中では今まで通りセル中心の値として使用します。

このような流れで全ての式に対しフィルタリング操作，ファーブル平均操作を施します．フィルタリング操作を施した質量保存式，運動量保存式，エネルギー保存式はそれぞれ次の通りです．

\begin{equation}
\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial }{\partial x_j}\left(\overline{\rho} \tilde{u_j}\right) = 0
\end{equation}

\begin{equation}
\frac{\partial \left(\bar{\rho} \tilde{u}_{j} \right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\bar{\rho} \tilde{u}_i \tilde{u}_j + \delta_{ij} \bar{p} \right) – \frac{\partial \breve{\sigma}_{ij}}{\partial x_j} = -\frac{\partial \tau_{ij}}{\partial x_j} + \frac{\partial }{\partial x_j} (\bar{\sigma}_{ij} – \breve{\sigma}_{ij})
\end{equation}

\begin{equation}
\frac{\partial \bar{e}}{\partial t} + \frac{\partial }{\partial x_j}\left[\left(\bar{e}+ \bar{p}\right) \tilde{u}_j \right] – \frac{\partial \breve{\sigma}_{ij} \tilde{u}_{i}}{\partial x_{j}} + \frac{\partial \breve{q}_{j}}{\partial x_{j}} =
-\frac{\partial }{\partial x_j}\left( C_{p} Q_{j} + J_{j} – D_{j} -(\overline{q}_{j}-\breve{q}_{j})\right)
\end{equation}

ただし、定義は下記の通りです。

\begin{eqnarray}
\bar{p} = \bar{\rho} R \tilde{T}
\end{eqnarray}

\begin{equation}
\breve{\sigma}_{ij} = \mu(\tilde{T}) \left(2 \tilde{S}_{ij} – \frac{2}{3} \delta_{ij} \tilde{S}_{kk} \right) \end{equation}

\begin{equation}
\tilde{S}_{ij} = \frac{1}{2} \left( \frac{\partial \tilde{u}_{j}}{\partial x_{i}} + \frac{\partial \tilde{u}_{i}}{\partial x_{j}}\right)
\end{equation}

\begin{equation}
\bar{e} = \bar{\rho} \tilde{E} = \frac{\bar{p}}{\gamma – 1} + \frac{1}{2} \bar{\rho} \tilde{u}_{i} \tilde{u}_{i} + \frac{\tau_{ij}}{2}
\end{equation}

\begin{equation}
\breve{q}_{j} = \lambda(\tilde{T}) \frac{\partial \tilde{T}}{\partial x_j}
\end{equation}

また，\(\tau_{ij}\)はSGS応力項，\(Q_{j}\)はSGS熱流束項，\(J_{j}\)はSGS乱流拡散項，\(D_{j}\)はSGS粘性拡散項です．

\begin{eqnarray}
\tau_{ij} = \bar{\rho} (\widetilde{u_{i} u_{j}} – \tilde{u_{i}}\tilde{u_{j}})
\end{eqnarray}

\begin{eqnarray}
C_{p} Q_{j} + J_{j} = (\overline{\rho u_{j} E} – \bar{\rho} \tilde{u}_{j} \tilde{E}) + (\overline{u_{j} p} – \tilde{u}_{j} \bar{p}) \end{eqnarray}

\begin{eqnarray}
J_{j} = \frac{1}{2} (\bar{\rho} \widetilde{u_{j} u_{i} u_{i}} – \bar{\rho} \tilde{u}_{j} \widetilde{u_{i} u_{i}})
\end{eqnarray}

\begin{eqnarray}
D_{j} = \overline{\sigma_{ij} u_{j}} – \breve{\sigma}_{ij} \tilde{u}_{j}
\end{eqnarray}

以上の式変形から、数値解析において必要なことを下記に記載します．

・運動量保存式の最後の項では，\(\bar{\sigma}_{ij} – \breve{\sigma}_{ij}=0\)を仮定する．
（粘性係数は温度の非線形関数のため，本来は0でないが影響が小さいため．）

・エネルギー保存式の最後の項では，\(\overline{q}_{j}-\breve{q}_{j}=0\)を仮定する．
（粘性係数は温度の非線形関数のため，本来は0でないが影響が小さいため．）

・\(\bar{\phi}\)や\(\tilde{\phi}\)は，コード上では\(\phi\)として扱い特別な操作をしない．セル中心で値を定義しているが，この値を代表としてセル内が一定の値であれば平均値も当然一緒なためである．

・\(D_{j}\)の影響は5%程度なため，無視する．

・3つのSGS項(\(\tau_{ij}, Q_{j}, J_{j}\))をモデル化し，計算できる形にする．

次回は3つのSGS項をモデル化していきたいと思います。

by hide

参考文献

Garnier, Eric and Sagaut, Pierre and Adams, Nikolaus,”Large Eddy Simulation for Compressible Flows”, Springer Netherlands, 2009

https://www.springer.com/gp/book/9789048128181

Pino Martín, M., Piomelli, U. & Candler, G. Subgrid-Scale Models for Compressible Large-Eddy Simulations . Theoret. Comput. Fluid Dynamics 13, 361–376 (2000). https://doi.org/10.1007/PL00020896