第13-3回　Smagorinsky モデル[julia]

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

前回はLESを考慮した支配方程式を導出しました。

支配方程式を解くうえで必要な、下記三つの値に関してモデル化していきます。

\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}

・SGS応力項(Smagorinskyモデル)

今回の解析では、一番有名なSmagorinskyモデルを使用します。

省略しますがこのモデルでは、GSからSGSへのエネルギーの輸送のみを考慮するモデルです。本来ならSGSからSGへのエネルギーの変換もありますが無視しています。

まずはBoussinesqの渦粘性近似を行います。

\begin{eqnarray}
\tau_{ij}-\frac{1}{3} \delta_{ij}\tau_{kk} = -2 \bar{\rho} \nu_{sgs}( \tilde{S}_{ij} -\frac{1}{3} \delta_{ij} \tilde{S}_{ij})
\end{eqnarray}

\begin{eqnarray}
S_{ij} = \frac{1}{2} \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) =
\left[\begin{array}{ccc}
\frac{\partial u}{\partial x} & \frac{1}{2} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & \frac{1}{2} \left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\ \frac{1}{2} \left(\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) & \frac{\partial v}{\partial y} & \frac{1}{2} \left(\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\ \frac{1}{2} \left(\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) & \frac{1}{2} \left(\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \right) & \frac{\partial w}{\partial z} \end{array}\right]
\end{eqnarray}

\(\tau\)も\(S_{ij}\)もテンソルであることに注意してください。

また，\(\tilde{S}_{ij}\)は GS のひずみ率テンソル(Strain-rate tensor)であり，\(\nu_{sgs}\)は乱流渦(動)粘性係数((kinematic) eddy viscosity)と呼ばれます。乱流渦粘性係数は場所と時間に依存して変わるため，毎回計算する必要があります．

SGS運動エネルギーの輸送方程式において，SGS運動エネルギーの生成と散逸がほぼ等しいと仮定します（局所平衡状態）

この仮定のもと、次元解析から乱流渦粘性係数\(\nu_{sgs}\)は次のように与えられます。

\begin{eqnarray}
\nu_{sgs} = (C_S \Delta)^2 |\tilde{S}| \
\Delta = \sqrt[3]{h_x h_y h_z} = \sqrt[3]{V} \
|\tilde{S}| = \sqrt{2 \tilde{S}_{ij} \tilde{S}_{ij}}
\end{eqnarray}

ここで，\(C_s=0.18\)です。また，テンソルの縮約規約から下記のように計算されます。

\begin{eqnarray}
S_{ij} S_{ij} &=& S_{11}S_{11} + S_{12}S_{12} + S_{13}S_{13} \\
&+& S_{21}S_{21} + S_{22}S_{22} + S_{23}S_{23} \\
&+& S_{31}S_{31} + S_{32}S_{32} + S_{33}S_{33}
\end{eqnarray}

解析時には、セル境界における\(\tau\)を計算する必要があります。

計算方法としては、粘性項の計算と同様に中心差分を用いて計算します。

・SGS熱流束項

SGS熱流束項をGS温度勾配の比と仮定し、導出します。

\begin{eqnarray}
\lambda_{sgs} = \frac{\mu_{sgs} C_{p}}{Pr_{sgs}}
\end{eqnarray}

\begin{eqnarray}
Q_{j} = -\frac{\bar{\rho} \nu_{sgs}}{Pr_{sgs}} \frac{\partial \tilde{T}}{\partial x_{j}}
\end{eqnarray}

ただし、下記の通りである。

\begin{eqnarray}
\nu_{sgs} = Smagorinsky \ model
\end{eqnarray}

\begin{eqnarray}
Pr_{sgs} = 0.6 \ \ \ (0.3 \ to \ 0.9)
\end{eqnarray}

・SGS乱流拡散項

無視しているケースも多いですが、今回はknightらのモデルに沿ってモデル化しています。

\begin{eqnarray}
J_{j} = \tau_{jk} \tilde{u}_k
\end{eqnarray}

以上で計算できる状態にモデル化できました。

ただし、一つ問題があります。

境界条件として、速度0の壁条件を課す場合、SGS応力項は0にならなければいけません。加えて壁近傍は渦粘性が小さくなるのでその効果を導入します。

このことを考慮するため、次回は壁法則について行います。

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

Knight, D., Zhou, G., Okong’o, N., and Shukla, V. (1998). Compressible large eddy simulation using unstructured grids. AIAA Paper
98-0535.