図1: ラバルノズル
・二次元ラバルノズルの計算(PASOLの効果)
Phoenicsの特徴の一つが，基本的に構造格子を採用しているという点である．その利点は03/08/21 Phoenicsとは？に示した通りだが，構造格子には苦手種目もある．意外なことに，解析的に流れ場が計算できる二次元ラバルノズルの計算が余り得意でない．今回は，このケースについて考察しよう．
・二次元ラバルノズル
ラバルノズルとは，図1のように途中にスロートを持つ内部流れで，圧縮性流体において超音速流れを得るために用いられる．一般に，ノズル出口における流速，圧力は圧縮波，膨張波が複雑に相互作用するため一様ではない．しかし，圧縮波と膨張波を相殺することで，ノズル出口において均一な平行流を得るための設計手法がいくつか知られている．
図1: ラバルノズル
参考文献：たとえば松尾一奏「圧縮性流体力学」，理工学社 pp. 256-258．
代表的な設計手法であり，解が解析的に求まるので入門者向けなのが，Foelschの手法である．ありがたいことに，詳しい資料がNASAにより無償提供されている．
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19930082268_1993082268.pdf
これを元に，Mathematicaを用いてノズル設計スクリプトを作成した．
ノズル設計スクリプト
このスクリプトが吐き出すのは，ノズルの形状を(x,y)のペアで出力したものにすぎない．これをPhoenicsが解釈可能な形状ファイルに変換する必要があるが，これが結構面倒くさい．私は，こういう作業はPerlを使い自動化することにしている．
ノズルデータ→形状データ変換スクリプト
このようにしてできたのが図2のような系である．
図2：Phoenicsで定義された二次元ラバルノズル
設計マッハ数は3，変曲点における広がり角(θ1)は0.2radと小さめに取る．実は，これはわざとで，θ1がνt/2に近いと，下のような現象は見られなくなる．
ノズルの壁面はスリップ境界で，流体はγ=1.4の非粘性流体とした．このとき，圧縮性流体の理論によれば，ノズル出口において流れに垂直な断面を取ったマッハ数分布は均一にならなくてはいけない．
・二次元ラバルノズルの計算
まず，上計算を，(20×100)の比較的粗いメッシュで計算してみた．メッシュの密度は系全体にわたり均一である．本来は，流速変化の多いところの密度を上げるべきであるが，比較のために敢えて愚鈍なメッシュの切り方を採用した．
case00のq1ファイル
############################################################################### ## Case00 ##################################################################### ############################################################## 26 August 2003 # ##Domain and Variables ## Title: Case00 M=3.0 base ## Mesh number x= 1 y= 20 z=100 ## Domain size 1×2.11×23.2 [mm] ## Flow: 0.714 - 0 [mmol/s] ## Tubulence model OFF ## Viscosity model Unknown ## Outlet pressure 400.0 [Pa] coeff. 1 ## ## Calculation options: ## EXPERT = ON ## DENPCO = ON ## WFIELD = P1 , V1 , W1 , ## SOLVE = SOLVE(P1 ,V1 ,W1 ,TEM1) ## STORE = STORE(DEN1,PRPS,MACH) ## ## Calculation convergence ## variable resref (res sum)/resref (res sum) ## P1 1.151E-09 1.579E+00 1.817E-09 ## V1 3.006E-08 6.370E+00 1.915E-07 ## W1 4.846E-07 9.615E-01 4.659E-07 ## TEM1 2.063E-04 2.389E-01 4.929E-05 ##
うーん，ノズル壁面にそって境界層(のようなもの)ができてしまった．構造格子を採用すると，原理的に"ちょうど壁面"の物理量を定義することができない．壁面上の流速は直近の格子点から推測するしかないわけだが，今回のケースのように，流速が零
(ノズルの中)から急激にM>1(気体)に変化するような場合，誤差が大きくなるのではないだろうか．
では，メッシュ密度を上げるとどうなるだろうか．メッシュ数を(100×180)に変更．密度は均一である．
############################################################################### ## Case01 ##################################################################### ############################################################## 26 August 2003 # ##Domain and Variables ## Title: Case01 M=3.0 High Mesh ## Mesh number x= 1 y=100 z=180 ## Domain size 1×2.11×23.2 [mm] ## Flow: 0.714 - 0 [mmol/s] ## Tubulence model OFF ## Viscosity model Unknown ## Outlet pressure 400.0 [Pa] coeff. 1 ## ## Calculation options: ## EXPERT = ON ## DENPCO = ON ## WFIELD = P1 , V1 , W1 , ## SOLVE = SOLVE(P1 ,V1 ,W1 ,TEM1) ## STORE = STORE(DEN1,PRPS,MACH) ## ## Calculation convergence ## variable resref (res sum)/resref (res sum) ## P1 2.854E-10 2.956E+01 8.438E-09 ## V1 8.596E-09 2.181E+01 1.874E-07 ## W1 1.244E-07 8.181E-01 1.017E-07 ## TEM1 5.322E-05 1.596E+00 8.496E-05 ##
相当に均一なマッハ数分布となった．ここから，このラバルノズルの流れを正確に計算するには(20×100)のメッシュでは全く不足であることがわかる．また，メッシュを充分細かく取れば，構造格子でも正確な計算が可能であることが示された．
・PASOLは使えるか?
Phoenicsでは，構造格子の欠点を補うために，「細分化格子」と「PASOL」という機能を提供している．これを使ってみよう．ただし，細分化格子はParallel版では使用できないので，我々としては頼りになるのはPASOLしかない．Case00の条件はそのまま，PASOLを適用して計算してみた．
############################################################################### ## Case02 ##################################################################### ############################################################## 26 August 2003 # ##Domain and Variables ## Title: Case02 M=3.0 base+PASOL ## Mesh number x= 1 y= 20 z=100 ## Domain size 1×2.11×23.2 [mm] ## Flow: 0.714 - 0 [mmol/s] ## Tubulence model OFF ## Viscosity model Unknown ## Outlet pressure 400.0 [Pa] coeff. 1 ## ## Calculation options: ## EXPERT = ON ## DENPCO = ON ## WFIELD = P1 , V1 , W1 , ## SOLVE = SOLVE(P1 ,V1 ,W1 ,TEM1) ## STORE = STORE(DEN1,PRPS,MACH) ## ## Calculation convergence ## variable resref (res sum)/resref (res sum) ## P1 1.145E-09 5.497E-01 6.293E-10 ## V1 4.587E-08 2.249E+01 1.032E-06 ## W1 4.690E-07 6.172E-01 2.895E-07 ## TEM1 2.118E-04 1.177E-01 2.493E-05 ##
駄目だった．Case00の計算が，そのままなめらかになっただけだった．つまり，「PASOLはメッシュ密度の不足を補うものではなく，メッシュが増えると結果が変わる系においてはPASOLを使うよりメッシュを増やせ」ということがわかる．
・現実の計算ではどうだ
今までの計算は，教科書では解析解が求まるものの，現実には存在しない流れである．だからこそ，ノズル表面でのマッハ数が3という答えが正解になる．では，実際に，現実の世界である「壁面ノンスリップ，粘性流体」という条件において，(20×100)と(100×180)の結果はどの程度異なるものだろうか．窒素ガスを流体として計算してみよう．今度は，壁面に境界層が発達し，マッハ数分布はなめらかに変化するので，構造格子の弱点が露わになることは無いのではないだろうか．
系のReynolds numberは，ノズル出口の流速，圧力で計算すると約340で，これはCOILの典型的なノズル下流のReと一致する．
case03のq1ファイル
case04のq1ファイル
############################################################################### ## Case03 ##################################################################### ############################################################## 28 August 2003 # ##Domain and Variables ## Title: Case04 <-03 M=3.0 nonslip, viscous ## Mesh number x= 1 y= 20 z=100 ## Domain size 1×2.11×23.2 [mm] ## Flow: 0.714 - 0 [mmol/s] ## Tubulence model OFF ## Viscosity model Sutherland's law ## Outlet pressure 400.0 [Pa] coeff. 1 ## ## Calculation options: ## EXPERT = ON ## DENPCO = ON ## WFIELD = P1 , V1 , W1 , ## SOLVE = SOLVE(P1 ,V1 ,W1 ,TEM1) ## STORE = STORE(DEN1,PRPS,MACH) ## ## Calculation convergence ## variable resref (res sum)/resref (res sum) ## P1 1.145E-09 1.553E+00 1.777E-09 ## V1 2.608E-08 6.761E+00 1.763E-07 ## W1 4.687E-07 9.233E-01 4.328E-07 ## TEM1 2.110E-04 2.790E-01 5.886E-05 ##
############################################################################### ## Case04 ##################################################################### ############################################################## 28 August 2003 # ##Domain and Variables ## Title: Case04 <-03 M=3.0 nonslip, viscous ## Mesh number x= 1 y=100 z=180 ## Domain size 1×2.11×23.2 [mm] ## Flow: 0.714 - 0 [mmol/s] ## Tubulence model OFF ## Viscosity model Sutherland's law ## Outlet pressure 400.0 [Pa] coeff. 1 ## ## Calculation options: ## EXPERT = ON ## DENPCO = ON ## WFIELD = P1 , V1 , W1 , ## SOLVE = SOLVE(P1 ,V1 ,W1 ,TEM1) ## STORE = STORE(DEN1,PRPS,MACH) ## ## Calculation convergence ## variable resref (res sum)/resref (res sum) ## P1 2.800E-10 3.147E+01 8.809E-09 ## V1 6.610E-09 3.012E+01 1.991E-07 ## W1 1.144E-07 1.015E+00 1.160E-07 ## TEM1 5.491E-05 9.766E-01 5.363E-05 ##
やっぱり，何となく合っていない．となると，流れに垂直な断面のマッハ数分布を正確に知るためには，相当細かくメッシュを刻まないとだめだなあ．
・まとめ
