１．はじめに
伝統的なマクロ衝撃波管についての研究は長い歴史があるが、近年の話題は、例えば マイクロ推進システムや医学システムにおける薬物伝達デバイスなどのような、マイクロ衝撃波管をより多くの産業分野で使用する試みである。
末端が閉じた衝撃波管は、低圧の非駆動ガスと高圧の駆動ガスは隔離されており、流れは隔膜が破裂することで発生する。
この破裂によって、低圧ガス中に衝撃波と衝撃波面が移動し、高圧ガスに膨張波が移動する、といった動きを結果として生じる。
マイクロとマクロの衝撃波管の主要な違いは、微細なサイズの流れの物理が必要になることである。特に、マイクロ衝撃波管では、低レイノルズ数における粘性の効果から衝撃波が減衰してしまうことがある。
一方、高クヌーセン数では、非連続体の影響から壁面近傍で流体がスリップし、これにより衝撃の大きさを増加させ、衝撃波の伝播を助けるように作動する。
以下は、マイクロ衝撃波管について、有名なCFDパッケージについて比較目的のために韓国安東国立大学の気体力学研究所で調査された幾つかの中の一つである。
我々CHAMは、図１で示すマイクロ衝撃波管の流れを計算するための設定方法などサポートをしました。
２．マクロ衝撃波管
マイクロ衝撃波管を計算する前段階として、我々CHAMは、解析解が既知である一般的なマクロ衝撃波管の一次元過渡計算を提案しました。
このシミュレーションの目的は、PHOENICS??が衝撃波、不連続面の接触、反射波等とらえられるかについての検証と性能評価である。
これらの現象は、マイクロ衝撃波管でも存在するが、マイクロ衝撃波管では大きな粘性効果を受け、高クヌーセン数による反射の効果も誘発される。
管の長さと経過時間は、２つの波は管の末端から反射する前に計算が終わるように選ばれなくてはならない。
シミュレーションは、0.1m2の断面積、長さ10mの管にたいし100等分割メッシュを使用した。
非定常計算は、6.1ミリ秒のシミュレーション時間を、100等分割した時間ステップで行われる。
初期条件は、速度は0m/sで個々の室では以下の設定を持つ。
駆動ガス： 圧力=1.0バール、温度= 348.391 K、および密度=1.0kg/m3
非駆動ガス： 圧力=0.1バール、温度= 278.13 K、および密度=0.125kg/m3
壁の摩擦は無視し、エネルギー方程式は理想気体の状態式によって静エンタルピーを計算する。
対流項に関し高次のVan Leer MUSCLスキームを使用し、時間項はデフォルトの一次オーダーのオイラースキームを使用する。
OpenFORM??を使用したCFDシミュレーションは、CHAMが最近リリースした「PH2OF」（PHOENICSデータからOpenFOAMデータを作成する）ユーザーインターフェイスを用いてOpenFOAM用のデータを自動的に作成した。
２つのデータはメッシュ、時間ステップ、対流項は同一なものである。
計算結果の精度に影響する様々な数値パラメータをいろいろ調整するしなくても、図２に示すように、PHOENICSとOpenFOAMおよび解析解の間で凡そ良い一致が見られた。
PHOENICSデフォルトのハイブリッドスキームでは、誤差のため誤った結果を引き起こすかもしれないので、Van Leer MUSCLスキームをこれらのシミュレーションに使用した。
図２に示すとおりOpenFOAMでは、計算結果の分布に対して多少、非物質的なアンダーシュートとオーバーシュートが発生した。
３．マイクロ衝撃波管
伝統的なマクロ衝撃波管では衝撃波および接触面は、本質的に管を通して一定速度で伝わる。
しかし、マイクロ衝撃波管では、減速するために、接触面を加速させる、はるかに厚い境界層が衝撃波の後ろに生成されることで衝撃波を減速させ、さらに、これら2つの波の間に不均一な流れが生じる。
図1で示すマイクロ衝撃波管の初期条件として、駆動、非駆動のチャンバー圧力は9.0と1.0気圧、温度300K、速度0m/sである。
このケースでは、非駆動ガスの圧力は、反射効果が無視できるほどに十分に高い。
この状況は、図３に示すとおり、安東大学により550μsec毎の静圧を測定する実験で調査された。図３は、以前示した図1におけるP1とP2と同義の、センサー位置をS1,S2の測定の結果と、Fluent??を使ってアンドン大学で行われたCFDの結果も含んでいる。
安東大学により、ここまでのすべての3つのコードを使用して、円筒座標の構造メッシュの非定常二次元軸対称シミュレーションが行われた。
この格子は、半径方向に160分割、軸方向に190分割（非駆動領域は80*1172、駆動領域は160*728）とした。
それぞれのシミュレーションは、500マイクロ秒間を、表1に示す均一な時間ステップを使用して計算された。
OpenFOAMでは、0.1より大きい時間ステップを使用した場合、収束させるのが難しかった。
それは、部分的にはPHOENICSがOpenFOAMより20%速く計算できたという理由になっている。
表1は、個々のCFDコードにより使われたエネルギー方程式の公式化と低レイノルズ乱流モデルを含む主要なコンピュータ処理の詳細もリストする。
追加として、OpenFOAMではJones -Launder低レイノルズ数k-eモデルを実施したが、これは結果にほとんど違いを生じさせなかった。
すべてのシミュレーションについて、サザーランド則による分子粘性、エネルギー方程式のため比熱1004のJ/kgK、および層流、乱流プラントル数はそれぞれ0.71と0.86を使用した。
PHOENICSとOpenFoamシミュレーションは、Intel Corei7 3.4GHzインテルコアi7 CPU、16GbyteRAMの環境で計算された。
図4は、PHOENICSおよびOpenFOAMで計算された、隔膜の破裂の後約0.1ミリ秒の密度、温度、圧力、および絶対速度のコンターを比較したものである。
分布は、非駆動領域に対して右に動いている正常な衝撃波、および駆動領域に対して左に動いている膨張波を示す。
衝撃波が、非駆動ガスに対して温度、圧力、および密度の急激な変化により引き起こされる激しい加速を与えるが見られる。
ここでは示されないが、衝撃波が、閉じられた終端の壁から後ろに反射される際、非駆動ガスの温度、圧力、および密度が計算結果よりも大きな増加を示すことも表している。
図 4に示すよう、OpenFOAMでは、非駆動領域の直径縮小部より下流でより鋭い波が発生するだけでなく、駆動領域に膨張する鋭い反射波を表している。
OpenFOAMのギザギザした波のキャプチャは、たぶんk-e SSTモデルで作られた、ずっと小さな乱流粘性のためである。
特に、2レイヤーk-eモデルは、衝撃波と接触断面を横切る大きな乱流粘性を作ってしまい、鈍った効果を引き起こしている。
Q1インプットファイルのINFORM機能によってPHOENICSのk-e SSTモデルを組み込むことがかなり有効的であろう。
PHOENICSの追加の試みは、高速で圧縮性の強い流れの膨張の効果のため、計算された乱流レベルを減らすために2レイヤーk-eモデルに対してサルカール圧縮率補正[3]を使うようにされたものである。この部分修正は結果に対しほとんどの変化をもたらさなかった。
図5と6は、非駆動領域の壁面に置かれた2つのセンサーポイントで、測定されたおよび計算結果の圧力履歴を比較したものである。
PHOENICS、OpenFOAM、Fluentの計算結果と実験的なデータの間で比較した。すべての3つのコードによる計算結果は、おおよそ実験と同じような傾向にある。
PHOENICSとFluentは、初期で、衝撃波の反射を表すS1とS2おける衝撃波到着時刻をそれぞれよくとらえられている。
しかし、OpenFOAMは、反射した衝撃波のために、より速い到着時間となった。
OpenFOAMとPHOENICSは、シミュレーションの遅い時間に、衝撃波管の非駆動領域内の圧力レベルが同程度の結果になっている。
現象時間内をとおして圧力に対しより良いレベルを予測できたのはFluentであることが明らかである。
たぶん二次精度の時間差分とずっと小さい時間ステップで計算したためと思われる。
しかし、Fluentで計算された圧力での大きい変動は、測定で存在しない。
それらの存在は 、線形風上法の無制限の性質によって引き起こされる非物理的なアンダーシュートとオーバーシュートからくるのかもしれない。
そのような非物理的はあらは、Lamnaouer [2]による大型の衝撃波管の彼のシミュレーションにおいても観察された。
それを彼はFluentではMUSCLスキームにみられる流束に制限がないことに起因すると考えた。
研究はさらに、3つのコードと実験の間の不一致の理由を調査することが必要である。
| 
