・[mm]系単位で計算する

・まえがき
　Phoenicsの新バージョン，3.6.0がこの春にリリースされた．CHAM社の説明では，変更はユーザーインターフェースの部分が主で，ソルバーは変わっていないと言う．しかし，3.5.1では計算できていたケースのうち，3.6.0で計算できなくなってしまったものが出た．

　CHAMに問い合わせたところ，

なお、先日の”パーソルで格子分割が増えると速度が計算されない件”は、
色々と試行錯誤したのですが、こちらで解決できなかったので
イギリスの開発チームに報告しました。

備考：
PHOENICSは単精度のプログラムなので、体系の小さな計算だと倍精度のプログラムに
比べると
必然的に誤差が大きくなってしまうでしょう。
今回のケースのような体系の計算をする場合は、mm単位系に換算して計算する方が安
全だと思います。
試しに、今回のモデルをmm単位系にしてテストして見ましたが、格子分割を増やして
も速度は問題なく計算されました。

との返事．悩んでいた系は，直径0.5mm程度のISOノズルだったので，なるほど，と思う．で，今回初めて[mm]系単位での計算に挑むことになった．

・単位系の変換
　CHAMのサポートページでやり方を見る．そして分かったのは，Phoenicsには別に単位系を変換するような特別な機能はないということ．流体力学には，「流れの相似則」という法則がある．レイノルズ数Reを

$\displaystyle Re = \frac{\rho v l}{\eta}$

(1)

[ $\rho$ ] 流体の密度
[] 代表長さ
[] 代表速度
[ $\eta$ ] 流体の粘性

と定義するとき，Reの同じ流れは全て同じに振る舞う，というものだ．流体の実験では普通に行われる手法(例えば，模型の船と水槽での試験)だが，これをPhoenicsの中でやるだけに過ぎない．いま，系の代表長さが0.1mmのオーダーなので，これを1,000倍に拡大することを考える．すると式(1)から，流体の密度を1/1,000とすれば，流れの速度は同じになる，というわけだ．結果的に，これは長さの単位を[mm]に，密度の単位を[mg/mm^3]に，すべて「ミリ」で考えたときに一致する．そのほかの単位の換算をまとめると以下のようになる．

物理量	Phoenics標準単位[SI]	[mm]系単位	備考
長さ	[m]	=10^3[mm]
速度	[m/s]	=1[mm/ms]
密度	[kg/m^3]	=10^-3[mg/mm^3]
圧力	[Pa]=[kg/m・s^2]	=10^-3[mg/mm・ms^2]
温度	[K]	=[K]	影響を受けない
動粘度	[m^2/s]	=10^3[mm^2/ms]

※層流，定常の系で必要な物理量のみ記す．

・例題
　今回取り上げたのは，ISO-9300標準ノズル．粘性の効果を顕著に出すために，敢えてDischarge coefficientの小さい領域で計算する．ISOノズルとDischarge coefficientについては「Discharge coefficientの計算」を参照のこと．

解析領域(y × z)	0.0026×0.007[m]
メッシュ	26×70
流体	窒素(理想気体近似)
オブジェクト境界	smooth wall friction
ノズル径	0.001[m]
質量流量	2.65e-6[kg/s]
Cd	約0.82

[m]系で素直に定義したq1ファイルはこちら．

　さて，ここからが本題．これを[mm]系に書き換えるにはどうするか．

系の大きさを1,000倍にする．これは，当たり前ですね．オブジェクトも1,000倍に拡大．見た目は全く同じだが大きさが1,000倍の系ができる．
入口境界条件の変更
- 入口境界は，このケースでは「密度」「流速」で定義されている．従って，密度を1/1000に変更．流速は換算係数1なのでそのまま．
出口境界条件の変更
- 出口境界は，「圧力」で与えられているので，これを1/1,000とする．
物性の変更
- 今回のケースの場合，変更は必要ない．密度については，定義式がρ=P/(GASCON・T)で，Pもρも[m]系の1/1000なので変更の必要がない．一方，粘性係数は理想気体近似の場合はSutherlandの式で定義する必要がある(この件についての大失敗は以前報告した)が，これもENUL=A・T^(1.5)/(B T ρ)で計算され，[mm]系のばあいENULは1,000倍，密度は1/1000なので定数Aの大きさはそのままでよい．
- 熱伝達は今回は計算しないので省略．CHAMホームページによれば，換算係数は1でいいそうだ．
初期条件の変更
- 初期条件は，普通の計算の場合domain fluidの「初期圧力」「初期温度」を定義する．初期圧力を1/1000にしておくのを忘れないように．
緩和係数の変更
- これは，以外と気がつかないのでは？Phoenicsの緩和係数は，圧力と密度についてはLinear relaxation coefficientで無次元量だが，その他の物理量はFalse time stepなので時間の次元を持つ．これを，1000倍にしないと計算が収束しない．こいつに気がつくまで苦労した．今回の計算の場合，V1，W1，T1の緩和係数を1,000倍にした．

このようにして出来た[mm]系のq1ファイルはこちら．両者を比較してみると良いでしょう．

・計算結果
　両者の計算結果を比較してみる．まずは，計算収束の状況．

[m]系

[mm]系

プローブ上の物理量はぴたりと一致した．収束の状況もうり二つ．圧力の単位が1,000倍違うことに注意．

続いて，速度コンター．
[m]系

[mm]系

全く区別がつきませんな．

最後に，中心軸上の物理量を比較してみる．

このグラフ，実線と点線で区別出来るようにしたのだが，重なって全く区別できない．

・まとめ

Phoenicsにおいて，長さの単位を[mm]にする方法について解説した．
Phoenicsには単位系を変更する機能はなく，単に「大きさ1,000倍，密度1/1,000」の流れ場を計算しているに過ぎない．
既にある[m]系の系を[mm]系に移行するためには，大きさ，密度，圧力だけでなく緩和係数も変えなくてはいけない．
[m]系と[mm]系の計算結果は計算精度の範囲内で一致．
ついでに，3.6.0でいままで出来なかったメッシュ数の計算を[mm]系でやらせてみたが，問題無くできた．