Turbulence Model / 乱流モデル

STAR CCM+で使う乱流モデルについて簡単にまとめてみました。ただ、多くのテキストは定義よりも、そのモデルのメリット・デメリットだけ述べているものが多いので、ここでは、理解して使うための定義にフォーカスしたまとめを簡単に記しています。

Preliminary

今更説明はいらないと思うけれど、大事だし、話を進めるうえでの前提条件となっているので一応… まず、Navier–Stokes equations は圧縮性、非圧縮性流体でそれぞれ次の通り書き表すことができる。

● 圧縮性流体 / Compressible flow

\frac{\partial ρ}{\partial t} + \frac{\partial}{\partial x_{i}} (ρ u_{i}) = 0

\frac{\partial}{\partial t} (ρ u_{i}) + \frac{\partial}{\partial x_{j}} (ρ u_{j} u_{i}) = - \frac{\partial p}{\partial x_{i}} + \frac{\partial τ_{i j}}{\partial x_{j}}

\frac{\partial}{\partial t} (ρ E) + \frac{\partial}{\partial x_{j}} (ρ u_{j} H) = \frac{\partial}{\partial x_{j}} (u_{i} τ_{i j}) + \frac{\partial}{\partial x_{j}} (k \frac{\partial T}{\partial x_{j}})

ただし、

ρ

と　

μ

は変数 (variables)。

● 非圧縮性流体 / Incompressible flow

\frac{\partial u_{i}}{\partial x_{i}} = 0

\frac{\partial u_{i}}{\partial t} + u_{j} \frac{\partial u_{i}}{\partial x_{j}} = - \frac{1}{ρ} \frac{\partial p}{\partial x_{i}} + \frac{1}{ρ} \frac{\partial τ_{i j}}{\partial x_{j}}

\frac{\partial T}{\partial t} + u_{j} \frac{\partial T}{\partial x_{j}} = k \nabla^{2} T

ρ

と　

μ

は定数 (constant)。

p

が運動量方程式 (momentum equation) だけにあることに気づいてね。

また、Boussinesq は乱流モデルを解くうえで、Laminar flow の場合、剪断応力は平均ひずみ速度 (mean strain rate) の線形モデルで近似する仮定が成立すると考えた (1877)。このモデルの比例係数が、いわゆる、eddy viscosity と呼ばれるもの。

● Navier–Stokes (or Laminar flow)

Viscous stress

τ_{i j} = 2 μ S_{i j} - \frac{2}{3} μ \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}

Strain sate

S_{i j} = \frac{1}{2} (\frac{\partial u_{k}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})

● RANS

Reynolds stress

τ_{i j}^{R} = 2 μ_{T} \overset{―}{S_{i j}} - \frac{2}{3} ρ K δ_{i j}

Reynolds averaged strain rate

\overset{―}{S_{i j}} = \frac{1}{2} (\frac{\partial \overset{―}{u_{i}}}{\partial x_{j}} + \frac{\partial \overset{―}{u_{j}}}{\partial x_{i}})

RANS k-ε model

RANS の二次方程式渦粘性モデルには、k-ε、k-ω などがあるけれど、おそらく、k-ε が現時点 (2020年) において、一般産業では最もよく使われている手法で、時間平均をベースに状態を計算するもの。

Dissipation rate

ρ \frac{\partial ϵ}{\partial t} + ρ v_{j} \frac{\partial ϵ}{\partial x_{j}} = C_{ϵ 1} \frac{ϵ}{k} τ_{i j} \frac{\partial v_{i}}{\partial x_{j}} - C_{ϵ 2} ρ \frac{ϵ^{2}}{k} + \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{T}}{σ_{ϵ}}) \frac{\partial ϵ}{\partial x_{j}}]

Turbulent kinetic energy

ρ \frac{\partial k}{\partial t} + ρ v_{i} \frac{\partial k}{\partial x_{j}} = τ_{i j} \frac{\partial v_{i}}{\partial x_{j}} - ρ ϵ + \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{T}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}]

ただし、k と

ϵ

は変数であり、

μ_{T} = C_{μ} ρ k^{2} / ϵ

NOTE

wall grid の要件を細かく指定することで、CFDの収束を良くすることができる。

Large Eddy Simulation (LES)

元々、LES の基本となる考え方は気象学のために、Smagorinsky によって提唱された (1963)。Smagorinsky-Lilly model は次の式で定義される。

v_{t} = (C_{S} Δ_{g})^{2} \sqrt{2 \overset{―}{S_{i j}} \overset{―}{S_{i j}}} = (C_{S} Δ_{g})^{2} | S |

こちらは時間平均ではなく、grid size よりも大きな渦だけ計算して、grid sizeよりも小さいもの (subgrid scales) はモデル化して計算する、filtering (と呼んでいるけれど、日本語的には、空間を平均化するイメージ) ベースの計算手法。
それから、LES を使う上で大事なのが、適切な dissipation rate (日本語だと散逸率かな？) を選択すること。

References

Henne, P. A. (1990) “Applied Computational Aerodynamics”, American Institute of Aeronautics and Astronautics

Scotti, A., Meneveau, C., Fatica, M. (1996) “Dynamic Smagorinsky model on anisotropic grids”, NASA Technical Reports, NASA-TM-111953.