乱流の渦構造可視化のためのQ値導出について
渦構造を可視化するのにQ値なるものがあるらしい．Q値って何？ってその人に聞いても教えてくれなかったので書き下し．流体をちょこっと触れていても乱流についてはそこまで触れてこなかったので，こういうところで差が出てきますね．知らないと圧力場とか渦度ベクトルの絶対値とかでやりきろうとするもんね．
流れとしては以下のような見つけ方です．
- 渦構造を特徴づけるのは各流体要素の速度勾配テンソル
- テンソル量を座標に依らないスカラーで取り扱いたい
- 速度勾配テンソルを対称テンソルと反対称テンソルに分けて変形速度テンソルと渦度テンソルに分解する
- それぞれの大きさを比較するとうまくテンソル不変量が出てくる
速度勾配テンソル
場の速度勾配を記述したテンソルを以下で定義できます，
\begin{align}
D= \begin{pmatrix}
\dfrac {\partial u}{\partial x} & \dfrac {\partial u}{\partial y} & \dfrac {\partial u}{\partial z}\\
\dfrac {\partial v}{\partial x} & \dfrac {\partial v}{\partial y} & \dfrac {\partial v}{\partial z}\\
\dfrac {\partial w}{\partial x} & \dfrac {\partial w}{\partial y} & \dfrac {\partial w}{\partial z}
\end{pmatrix}\end{align}
ここでは座標で
はその座標系での速度です．なぜこれを定義するかといえば，各点(厳密に言えば流体要素)の渦構造を特徴づけるのは速度の絶対値ではないのでその勾配をとってみようということではないかと思います．多分3階4階のテンソルが有用ならそっちにも言っていると思いますが，多分2階テンソル量である程度有用なのでしょう．ちなみに速度の絶対値は相似則を考えればわかるように渦構造とはリンクしていないでしょう．最初は渦度じゃダメなの？と思いましたが渦度は単純せん断場でも値を持ちますし，ベクトル量なので扱いづらいし座標系によって変わっちゃうから？(変わるかどうか調べていません)と思います．
テンソル不変量
テンソルには座標系に依らないスカラー量が定義できます．そのうち第二不変量というのを今回使うのでそれをとすると
の成分表示で次のように表せます．
\begin{eqnarray}
I_2=\begin{vmatrix}
D_{22} & D_{32}\\
D_{23} & D_{33}
\end{vmatrix}+\begin{vmatrix}
D_{11} & D_{21}\\
D_{12} & D_{22}
\end{vmatrix}+\begin{vmatrix}
D_{11} & D_{31} \\
D_{13} & D_{33}
\end{vmatrix}
\end{eqnarray}
それらの導出は以下に超詳しく書いてあるのでこちらをご覧ください．
これらの量に物理的な意味を持たせられたらスカラー値で取り扱えてしかも座標系に依らないものが手に入るよね，ということでしょう．
対称テンソルと反対称テンソル
まずはこの速度勾配テンソルを対称テンソルと反対称テンソルに分解してみます.
\begin{eqnarray}
D_{ij}=\dfrac {1}{2} \left( D_{ij}+D_{ji} \right) +\dfrac {1}{2} \left( D_{ij}-D_{ji} \right)
\end{eqnarray}
ここで右辺の各項は流体要素の変形成分である変形速度テンソルと回転成分である渦度テンソル
にあたります。
\begin{eqnarray}
S_{ij}=\dfrac{1}{2}\left( D_{ij}+D_{ji} \right)\\
\Omega_{ij} = \dfrac{1}{2}\left( D_{ij}-D_{ji} \right)
\end{eqnarray}
というわけで速度勾配テンソルから変形とか渦とか物理的意味のあるものが出てきました。この量とテンソル不変量が結びつけばQ値の物理的意味がわかりやすいですね。
テンソル不変量で表してみる
さて、少し天下り的ですが、各テンソル同士のスカラー積(複内積というそうです)の差，なる量を計算してみます。これは物理的には速度場の回転成分と変形成分の差のようなものでしょう。
\begin{align}
\Omega ^{2}_{ij}-S^{2}_{ij} &=\left( \Omega _{ij}+S_{ij} \right) \left( \Omega _{ij}-S_{ij} \right)\\
&=-D_{ij}D_{ji}\\
&=-(D_{11}D_{11}+D_{12}D_{21}+D_{13}D_{31}+D_{21}D_{12}+D_{22}D_{22}+D_{23}D_{32}+D_{31}D_{13}+D_{32}D_{23}+D_{33}D_{33})\\
&=-\left( D_{11}+D_{22}+D_{33} \right) ^{2}+2 \left( D_{11}D_{22}+D_{22}D_{33}+D_{33}D_{11} \right)
-2 \left( D_{12}D_{21}+D_{13}D_{31}+D_{23}D_{32} \right)\\
&=-2 \left( D_{12}D_{21}+D_{13}D_{31}+D_{23}D_{32}-D_{11}D_{22}-D_{22}D_{33}-D_{33}D_{11} \right)\\
&=2\left( \begin{vmatrix}
D_{22} & D_{32}\\
D_{23} & D_{33}
\end{vmatrix}+\begin{vmatrix}
D_{11} & D_{21}\\
D_{12} & D_{22}
\end{vmatrix}+\begin{vmatrix}
D_{11} & D_{31} \\
D_{13} & D_{33}
\end{vmatrix} \right) \\
&= 2I_2
\end{align}
とこれは不変量で表せましたね！
あとは係数の問題なのですがQ値としては以下の定義がされるみたいです．
\begin{eqnarray}
Q=\dfrac{1}{2}I_2=\dfrac{1}{2}\left( \Omega^2 - S^2 \right)
\end{eqnarray}
めでたしめでたし。
複 雑 内 部 流 れ 場 に お け る渦 構 造 の 同 定
https://www.jstage.jst.go.jp/article/jvs1990/22/1Supplement/22_1Supplement_107/_pdf
http://solid4.mech.okayama-u.ac.jp/テンソル.pdf
質問コメントなどありましたらご連絡いただけるとうれしいです．
しかしテンソルなんてあまり気にしてなかったけどちゃんとやってみるものだなあ．．．
メッシを切ってみる
というわけでこの前作ったジオメトリのメッシュを切ってみましょう。
メッシュを切るのはあまり困難がないはずなのですが、実は以前作成したジオメトリが悪かったらしく、なぜか流体と便器のメッシュが連続していませんでした。なので何度かトライ&エラーでここまで実は来ています。本当はこういうところも残した方がいいと思うんだけど、まずは進んで詳しくはまたそのうちと言うことにしましょう。Fluentは不連続なメッシュでも境界をきちんと設定すれば計算できるはずなんだけど、基本的にジオメトリーが連続しているならばメッシュも連続させたいよねということで、ジオメトリレベルで工夫してました。
メッシュは一体どの程度切ればいいかよくわからないけど、デフォルトでダメだったらもう少し細くすると言う感じで行ってみましょう。マシンパワーがないしね、できるとこから。
配色の関係でなんか便器が消えている感じに見える。。。
次にめんどくさいのは名前選択ですね。僕は割と後から境界条件を増やしたくなる性格なのでその度にまた名前選択をやり直すのは非常にめんどくさいです。まぁ今回はジオメトリーがそんなに複雑じゃないからまずこんな感じでやってみて後で必要だったら増やしてくっていうか感じでいいんじゃないですかね。
というかそもそも境界条件はどうしようかと言うのは決めてなかったなあ、次はそれらしい条件を決めなければいけない。
Design Modelerでジオメトリを作成する
というわけで解析対象を作るべくまずはデザインモデラーを使ってお絵描きです．実は数値計算するにあたってこの作業が1番めんどくさいと僕は思います．ここら辺の作業を人工知能( )が担ってくれればいいのになぁと思います．とまあそんなこと言っていても仕方ないのでまずは適当にお絵描きをしてみます．
多分計算資源の関係で3-Dの計算はきついと思うのでまずは2次元から始めてみましょう．普段3-Dの描画にはスペースクレームを使うんですが、いかんせん2次元のスペースクレームの使い方がよくわからないのでまだデザインモデラーからです．
で適当にスケッチからサーフェスを作成していくと、よくわからないエラーが出てしまう。これが出る根本原因とかを本当は調べた方がいいと思うんだけどめんどくさいのでいつもこういう時はラインボディーを作ってから辺からサーフェスで何とか作ってしまいます。
こんな感じに。ただこうするとなんとなくちょっとこのツリーが汚くなるんだよなぁ。そして後々へサーフェスの被りがあるからそれをブリアン演算で引き算したりしなくてはいけない。
とまぁ色々と面倒なことが起きてしまったけれどもひとまずジオメトリーは完成しました。
さて今回の解析対象は見ての通りうちのトイレです．最近寒いのでトイレと脱衣所用にオイルヒーターを買ったんですがその効果を見るためにちょっとした計算をしてみましょう．部屋に対してトイレの便器がとてつもなく大きいのは2次元だからです、多分。
