工学系大学院生のブログ

2人の大学院生による雑記ブログ

第13-4回 乱流粘性(y+,壁関数)[julia]


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



前回はSmagorinskyモデルを用いて、SGS応力項等をモデル化しました。


しかし、速度0の境界条件すなわち滑りなし壁条件を課した場合は、SGS応力項等は0にならなければいけません。加えて壁近傍の渦粘性を小さくしなければいけません。



そこで、乱流渦粘性係数に減衰関数\(f(y^{+})\)を乗じます。

\begin{eqnarray}
\nu_{sgs} = (C_S \Delta f(y^{+}) )^2 |\tilde{S}|
\end{eqnarray}


この減衰関数は、壁面近傍では\(f(y^{+})=0\)に漸近し、離れれば\(f(y^{+})=1\)に漸近する関数になります。


具体的には、下記のVan Driest型減衰関数を適用します。

\begin{eqnarray}
f(y^+) = 1-exp(-y^+/26)
\end{eqnarray}





さて、これで境界条件の問題は解決しましたが、\(y^{+}\)をどのように求めればいいのでしょうか。


\(y^{+}\)は壁からの無次元座標であり、壁面近傍の無次元化されたもろもろの値は下記の通りです。


\begin{eqnarray}
y^+ = \frac{u_{\tau} y}{\nu} \\
u^+ = \frac{u}{u_{\tau}} \\
u_{\tau} = \sqrt{\frac{\tau_w}{\rho}} \\
\tau_w = \mu\frac{\partial u}{\partial y}|_{wall}
\end{eqnarray}


壁近傍の平均速度分布は,壁座標(\(y^+\))と無次元化速度(\(u^+\))による関数として表されます.

\begin{eqnarray}
F(y^+, u^+) = 0
\end{eqnarray}




仮に壁面近傍の速度,壁からの距離,動粘性係数が与えられたら,\(F(y^+, u^+) = 0\)から非線形方程式が構成でき,\(u_{\tau}\)が求まります.


\(u_{\tau}\)がもとまれば定義から,それぞれの値がわかり,乱流渦粘性係数を算出できます.


壁法則の1つとして,Spalding則を紹介します.

\begin{eqnarray}
F(y^+, u^+) = u^+ – y^+ – e^{-k B} \left(e^{-k u^+} – 1 – (k u^+)-\frac{1}{2}(k u^+)^2-\frac{1}{6}(k u^+)^3 \right) \\
k = 0.4 \\
B = 5.5
\end{eqnarray}




この壁法則は、粘性低層も含んだ近似曲線であるため壁面第一格子点がどの層に合っても問題ありません。

(この辺りの説明は省くので、知らないという方は、乱流境界層、壁法則等で調べてみてください。)



解き方としては、Spalding則に対してニュートン法を適用します.



まず,\(u_{\tau}\)の式に変形します。

\begin{eqnarray}
F(u_{\tau}) = \frac{u}{u_{\tau}} – \frac{u_{\tau} y}{\nu} – exp(-k B) \left(exp(-k \frac{u}{u_{\tau}}) – 1 – (k \frac{u}{u_{\tau}})-\frac{1}{2}(k \frac{u}{u_{\tau}})^2-\frac{1}{6}(k \frac{u}{u_{\tau}})^3 \right)
\end{eqnarray}

ここで,\(u, y, \nu, k, B\)は既知です.


上式の微分は下記の通りです.

\begin{eqnarray}
F'(u_{\tau}) = -\frac{u}{u_{\tau}^2} – \frac{y}{\nu} + exp(-k B) \left(k \frac{u}{u_{\tau}^2} \right) \left(exp(-k \frac{u}{u_{\tau}}) – 1 – (k \frac{u}{u_{\tau}})-\frac{1}{2}(k \frac{u}{u_{\tau}})^2 \right)
\end{eqnarray}

これらの結果から,\(u_{\tau}=0.1\)等の適当な値から始め,下記の式が収束するまで計算します.

\begin{eqnarray}
u_{\tau, n+1} = u_{\tau, n} – \frac{F(u_{\tau, n})}{F'(u_{\tau, n})}
\end{eqnarray}

また収束判定として,無次元化したノルム1が\(10^{-4}\)以下となった時としました.

\begin{eqnarray}
\frac{u_{\tau, n+1} – u_{\tau, n}}{u_{\tau, n+1}} < 10^{-4}
\end{eqnarray}




以上で乱流項は終わりです。


次回は、これまでのまとめとしてコードに落としこみたいと思います。


by hide

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です