・圧縮性流体の計算(レーザー加工ノズルまわりの流れ)


　典型的な圧縮性流体の例題として，2002年度に行ったレーザー加工ノズルまわりの流れを取り上げる．といっても，このページはむしろ私自身の備忘録的な目的で作ったので，多少ずぼらな方法でモデルを作成している．まあ，それも含めてご参考ということで．

・SEMモデルから圧縮性流体モデルへ

　自由表面モデルへの挑戦で行ったレーザー加工ノズルまわりの流れを計算して，これを2002年度の圧縮性流体のモデルと比較しようとして，当時の計算が今回のSEMモデルと直接比較しにくいことに気がついた．理由は，系の定義サイズが違うこと，メッシュの切り方，カーフ端点のノズルとの相対位置等々が違うことなどだが，これを機会に，いちど圧縮性流体モデルの計算を一から構築し，SEMモデルと比較可能な計算をやり直そうと考えた．

　ここで大切なのは，計算モデルを変えるときは，ケースをはじめから作り直すことだ．VR Viewerによるモデル作成は一見簡単だが，中ではPhoenicsソルバに与える様々なパラメータを自動生成している．残念ながら，VR Viewerの制限で，モデルを変えると必要なくなるパラメータは，自動的に消えてくれない．たとえば，モデルを自由表面onからoffに変えても，変数VFOLは残ったままになる．他にも，リセットされないプロパティは幾つかあるのだが，とにかく，圧縮性流体の計算をやりたければ，白紙の状態からもう一度ケースを構築する必要がある．

・計算領域，形状ファイル，メッシュをコピーする
　しかし，今回，やりたい計算は，SEMモデルと同じ大きさの系，同じオブジェクトなので，これを流用しない手はない．まず，新しいケースを作るために，白紙のアイコンをクリック．ケースの名前だけを決めて，いったんVR Editorを終了させる．


　続いて，SEMモデルのq1ファイルを開き，以下のブロックを今作ったq1ファイルに上書きする．

  1. ブロック1: 系の大きさとメッシュ
    Groups 3, 4, 5 Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
    RSET(M,80,70,95,1.000000E-05)
  2. ブロック2: メッシュ，領域
    > DOM, SIZE, 3.000000E-02, 1.500000E-02, 2.000000E-02
    > DOM, MONIT, 1.514380E-02, 7.500000E-06, 1.001850E-02
    > DOM, SCALE, 1.000000E+00, 1.000000E+00, 1.000000E+00
    > DOM, SNAPSIZE, 1.000000E-02
    > GRID, RSET_X_1, 28,-1.500000E+00
    > GRID, RSET_X_2, 8, 1.000000E+00
    > GRID, RSET_X_3, 8, 1.000000E+00
    > GRID, RSET_X_4, 8, 1.000000E+00
    > GRID, RSET_X_5, 28, 1.500000E+00
    > GRID, RSET_Y_1, 10, 1.000000E+00
    > GRID, RSET_Y_2, 20, 1.000000E+00
    > GRID, RSET_Y_3, 40, 1.500000E+00
    > GRID, RSET_Z_1, 1, 1.000000E+00
    > GRID, RSET_Z_2, 17,-1.500000E+00
    > GRID, RSET_Z_3, 40,-1.300000E+00
    > GRID, RSET_Z_4, 15, 1.000000E+00
    > GRID, RSET_Z_5, 21, 1.500000E+00
    > GRID, RSET_Z_6, 1, 1.000000E+00
  3. ブロック3: オブジェクト
    > DOM, SIZE, 3.000000E-02, 1.500000E-02, 2.000000E-02
    > DOM, MONIT, 3.716320E-04, 7.500000E-06, 1.000000E-04
    > DOM, SCALE, 1.000000E+00, 1.000000E+00, 1.000000E+00
    > DOM, SNAPSIZE, 1.000000E-02
    > GRID, RSET_X_1, 28,-1.500000E+00
    > GRID, RSET_X_2, 8, 1.000000E+00
    > GRID, RSET_X_3, 8, 1.000000E+00
    > GRID, RSET_X_4, 8, 1.000000E+00
    > GRID, RSET_X_5, 28, 1.500000E+00
    > GRID, RSET_Y_1, 10, 1.000000E+00
    > GRID, RSET_Y_2, 20, 1.000000E+00
    > GRID, RSET_Y_3, 40, 1.500000E+00
    > GRID, RSET_Z_1, 1, 1.000000E+00
    > GRID, RSET_Z_2, 17,-1.500000E+00
    > GRID, RSET_Z_3, 40,-1.300000E+00
    > GRID, RSET_Z_4, 15, 1.000000E+00
    > GRID, RSET_Z_5, 21, 1.500000E+00
    > GRID, RSET_Z_6, 1, 1.000000E+00

    > OBJ, NAME, SLAB1
    > OBJ, POSITION, 0.000000E+00, 1.500000E-04, 1.000000E-02
    > OBJ, SIZE, 3.000000E-02, 1.485000E-02, 1.000000E-03
    > OBJ, CLIPART, cube14
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, BLOCKAGE
    > OBJ, COLOR-MODE, 2
    > OBJ, COLOR-VAL, 1
    > OBJ, MATERIAL, 100

    > OBJ, NAME, SLAB2
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 1.000000E-02
    > OBJ, SIZE, 1.490000E-02, 1.500000E-04, 1.000000E-03
    > OBJ, CLIPART, cube14
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, BLOCKAGE
    > OBJ, COLOR-MODE, 2
    > OBJ, COLOR-VAL, 1
    > OBJ, MATERIAL, 100

    > OBJ, NAME, LIQUID
    > OBJ, POSITION, 1.490000E-02, 0.000000E+00, 1.000000E-02
    > OBJ, SIZE, 3.000000E-04, 1.500000E-04, 1.000000E-03
    > OBJ, CLIPART, cube14
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, BLOCKAGE
    > OBJ, MATERIAL, 100

    > OBJ, NAME, BOTTOM
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 3.000000E-02, 1.500000E-02, 0.000000E+00
    > OBJ, CLIPART, cube12t
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, OUTLET
    > OBJ, PRESSURE, 0.000000E+00
    > OBJ, TEMPERATURE, SAME
    > OBJ, COEFFICIENT, 1.000000E+03
    > OBJ, VELOCITY, SAME , SAME , SAME
    > OBJ, TURBULENCE, SAME , SAME

    > OBJ, NAME, TOP
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 2.000000E-02
    > OBJ, SIZE, 3.000000E-02, 1.500000E-02, 0.000000E+00
    > OBJ, CLIPART, cube12t
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, OUTLET
    > OBJ, PRESSURE, 0.000000E+00
    > OBJ, TEMPERATURE, SAME
    > OBJ, COEFFICIENT, 1.000000E+03
    > OBJ, VELOCITY, SAME , SAME , SAME
    > OBJ, TURBULENCE, SAME , SAME

    > OBJ, NAME, WEST
    > OBJ, POSITION, 0.000000E+00, 1.500000E-02, 0.000000E+00
    > OBJ, SIZE, 3.000000E-02, 0.000000E+00, 2.000000E-02
    > OBJ, CLIPART, cube12t
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, OUTLET
    > OBJ, PRESSURE, 0.000000E+00
    > OBJ, TEMPERATURE, SAME
    > OBJ, COEFFICIENT, 1.000000E+03
    > OBJ, VELOCITY, SAME , SAME , SAME
    > OBJ, TURBULENCE, SAME , SAME

    > OBJ, NAME, SOUTH
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 0.000000E+00, 1.500000E-02, 2.000000E-02
    > OBJ, CLIPART, cube12t
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, OUTLET
    > OBJ, PRESSURE, 0.000000E+00
    > OBJ, TEMPERATURE, SAME
    > OBJ, COEFFICIENT, 1.000000E+03
    > OBJ, VELOCITY, SAME , SAME , SAME
    > OBJ, TURBULENCE, SAME , SAME

    > OBJ, NAME, NORTH
    > OBJ, POSITION, 3.000000E-02, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 0.000000E+00, 1.500000E-02, 2.000000E-02
    > OBJ, CLIPART, cube12t
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, OUTLET
    > OBJ, PRESSURE, 0.000000E+00
    > OBJ, TEMPERATURE, SAME
    > OBJ, COEFFICIENT, 1.000000E+03
    > OBJ, VELOCITY, SAME , SAME , SAME
    > OBJ, TURBULENCE, SAME , SAME

    > OBJ, NAME, HI_NZL
    > OBJ, POSITION, 3.000000E-03, 0.000000E+00, 1.250000E-02
    > OBJ, SIZE, 2.400000E-02, 1.200000E-02, 7.500000E-03
    > OBJ, CLIPART, upnzl
    > OBJ, ROTATION24, 15
    > OBJ, GRID, 2
    > OBJ, TYPE, BLOCKAGE
    > OBJ, MATERIAL, 100

    > OBJ, NAME, HI_PLG
    > OBJ, POSITION, 1.270000E-02, 0.000000E+00, 1.980000E-02
    > OBJ, SIZE, 4.600000E-03, 2.300000E-03, 2.000000E-04
    > OBJ, CLIPART, hfcyl
    > OBJ, ROTATION24, 1
    > OBJ, GRID, 2
    > OBJ, TYPE, BLOCKAGE
    > OBJ, MATERIAL, 100

    > OBJ, NAME, HI_IN
    > OBJ, POSITION, 1.270000E-02, 0.000000E+00, 1.980000E-02
    > OBJ, SIZE, 4.600000E-03, 2.300000E-03, 0.000000E+00
    > OBJ, CLIPART, hfcyl
    > OBJ, ROTATION24, 1
    > OBJ, GRID, 2
    > OBJ, TYPE, INLET
    > OBJ, DENSITY, 1.821000E+00
    > OBJ, VOLUFLOW, 4.000000E-04
    > OBJ, TEMPERATURE, 2.980000E+02
    > OBJ, INLET_C1, 1.000000E+00
    > OBJ, TURB-INTENS, 5.000000E+00

    > OBJ, NAME, LO_NZL
    > OBJ, POSITION, 1.660000E-02, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 1.000000E-02, 1.000000E-02, 1.000000E-02
    > OBJ, CLIPART, lonzl
    > OBJ, ROTATION24, 2
    > OBJ, GRID, 2
    > OBJ, TYPE, BLOCKAGE
    > OBJ, MATERIAL, 100

    > OBJ, NAME, LO_PLG
    > OBJ, POSITION, 1.660000E-02, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 1.000000E-02, 1.000000E-02, 2.000000E-04
    > OBJ, CLIPART, loin
    > OBJ, ROTATION24, 2
    > OBJ, GRID, 2
    > OBJ, TYPE, BLOCKAGE
    > OBJ, MATERIAL, 100

    > OBJ, NAME, LO_IN
    > OBJ, POSITION, 1.660000E-02, 0.000000E+00, 2.000000E-04
    > OBJ, SIZE, 1.000000E-02, 1.000000E-02, 0.000000E+00
    > OBJ, CLIPART, loin
    > OBJ, ROTATION24, 2
    > OBJ, GRID, 2
    > OBJ, TYPE, INLET
    > OBJ, DENSITY, 2.000000E+00
    > OBJ, VOLUFLOW, 8.000000E-04
    > OBJ, TEMPERATURE, 2.980000E+02
    > OBJ, INLET_C2, 1.000000E+00
    > OBJ, TURB-INTENS, 5.000000E+00
    > OBJ, OBJECT-SIDE, 1.000000E+00

    > OBJ, NAME, WIRE_YZ
    > OBJ, POSITION, 1.400000E-02, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 2.000000E-03, 1.500000E-02, 2.000000E-02
    > OBJ, CLIPART, wirexyz
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, NULL

    > OBJ, NAME, WIRE_XZ
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 3.000000E-02, 1.000000E-03, 2.000000E-02
    > OBJ, CLIPART, wirexyz
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, NULL

    > OBJ, NAME, WIRE_XY
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 4.999999E-03
    > OBJ, SIZE, 3.000000E-02, 1.500000E-02, 5.000000E-03
    > OBJ, CLIPART, wirexyz
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, NULL

    > OBJ, NAME, WIRE_L
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 0.000000E+00
    > OBJ, SIZE, 3.000000E-02, 1.500000E-02, 2.000000E-04
    > OBJ, CLIPART, wirexyz
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, NULL

    > OBJ, NAME, WIRE_T
    > OBJ, POSITION, 0.000000E+00, 0.000000E+00, 1.980000E-02
    > OBJ, SIZE, 3.000000E-02, 1.500000E-02, 2.000000E-04
    > OBJ, CLIPART, wirexyz
    > OBJ, ROTATION24, 1
    > OBJ, TYPE, NULL

　これで，いったんq1ファイルをセーブする．続いて，オブジェクト定義で使われる，デフォルトでない形状ファイル，上ではLOIN.DAT，LONZL.DAT，HINZL.DATをq1と同じフォルダにコピー．これを忘れると，オブジェクトが豆腐になってしまう．

・計算条件の作成

　再びVR Editorを開けると，メッシュとオブジェクトがしたの図の様に再現されているはずである．

　まずは，Modelsの設定．圧縮性流体，高レイノルズ数流れの場合は下のようにする．エネルギー方程式は，静温度で解いた方が，圧縮性流体の場合には正確(CHAM社フォーラムより)．乱流モデルは必須であるが，どの乱流モデルを使うとより実際に近い流れ場を与えるかというのは，数値流体の最難問の一つ．ここは，素直にk-εモデルを使いましょう．

　続いて，このパネルのSolution cotrolサブパネルで，変数C1，C2を定義する．これは，質量を持たない物性量で，上ガス，下ガスをマークするために使われる．PhoenicsのSingle equation policyに従えば，任意の変数についてsolveの属性を与え，「連続の式」のみをその物性と定義すれば，その変数がマーカーになることは自明であろう．方法は，SOLVEのテキストボックスに「C1」と入力し，「Apply」ボタンで確定する．


上手くいくと，変数リストにC1とC2が現れる．


　Tips．高圧縮流体の場合，U，V，WはWhole fieldで解いた方が断然収束性がよい．

　続いて，Propertiesの設定を行う．ガスは，理想気体近似の窒素としてモデル化する．


　ここで重要なポイント．「Edit domain material properties」で，粘性モデルを必ず変更する必要がある．VR Viewerが作るモデルはデフォルトでは「動粘性が定数」で，流体を圧縮性としてもこれが変更されない．動粘性はガスの圧力に反比例するので，圧縮性ガスの場合は全く近似が成立しない．このヘタレな仕様のおかげで，導入当初は計算が合わずに随分苦労した．粘性モデルはいくつも用意されているが，予め物性値が呼び出されるのはSutherlandモデルだけ．その他の場合は自分で調べる必要があるので，今回はこれを利用する．


　続いて，PropertiesパネルからPage Dnで下がり，PILコマンドパネルを開く．ここで，C1とC2のプラントル数をセットする．いま，C1とC2は窒素ガスの広がりを表すマーカー変数としたいので，窒素と同じ拡散定数を与えてやれば拡散で広がる効果も模擬できる．窒素のプラントル数は0.716なので，以下のように入力．

　Numericsの設定．ループ回数は1,500とした．ここは，試してみて充分収束する回数を指定すればよい．

　Relaxationは，今回はこんな感じ．これも，経験と勘がモノを言う世界． 

 RELAX(P1  ,LINRLX, 1.000000E+00)
 RELAX(U1  ,FALSDT, 1.000000E-03)
 RELAX(V1  ,FALSDT, 1.000000E-03)
 RELAX(W1  ,FALSDT, 1.000000E-03)
 RELAX(KE  ,LINRLX, 5.000000E-01)
 RELAX(EP  ,LINRLX, 5.000000E-01)
 RELAX(TEM1,FALSDT, 1.000000E-03)

　

・境界条件の作成

　上ノズルの境界条件は下のようにする．ガスの密度をUser-setとして，ノズル内部のガス密度を与える．体積流量は，実験に合うように決めた．Setting scalerで，上ガスとともに流れ込むC1を1，C2を0とする．これで，上ガスをC1でマーキングすることになる．温度は，Reference temperatureからの相対値であることに注意．いまは，Reference temperatureは特にいじっていないので273Kである．従って，25.0指定は常温ということ．

　下ガスも同様．ただし，Setting scalarはC1=0，C2=1とする．

　流出境界．実は，この設定が重要である．詳しくはこちら．流出境界は，主に圧力で規定され，他の条件が計算結果に与える影響は比較的少ない．圧力の条件は当然0として，今回はその他の条件は以下の方針で決定．


全ての流出境界で同じ設定が出来たかをVR Viewerで確認するのは大変なので，q1ファイルを直接見る方がよろしい．設定も，copy&pasteでらくらくと出来るのは言うまでもない．

> OBJ, NAME, NORTH
> OBJ, POSITION, 3.000000E-02, 0.000000E+00, 0.000000E+00
> OBJ, SIZE, 0.000000E+00, 1.500000E-02, 2.000000E-02
> OBJ, CLIPART, cube12t
> OBJ, ROTATION24, 1
> OBJ, TYPE, OUTLET
> OBJ, PRESSURE, 0.000000E+00
> OBJ, TEMPERATURE, 2.500000E+01
> OBJ, COEFFICIENT, 1.000000E+00
> OBJ, TURBULENCE, SAME , SAME

これで，全てできあがり．


・VR EditorまたはEarthにバグがあることを発見した

　計算させてみたところ，問題があることが発覚した．slab1とslab2の間に隙間が空いてしまうのだ．VR Viewerでは，二つのオブジェクトはくっついているように見える．しかし，計算結果をVR Editorでみてみると，ちょうどメッシュ一つ分の空間があいていることがわかる．

どうやら，左側のslab2の境界を判断するときに，四捨五入で切り捨てられてしまったらしい．実は，浮動小数点表現の0.15は，2進数では循環小数で表される．したがって，四捨五入すると，最後の桁が切り捨てられ，0.149999・・・・になる可能性が充分ある．これをバグと言ってしまうのもかわいそうだが，q1ファイル上では二つのオブジェクトはくっついているはずなので，やはりバグとしか言いようがない．原因がVR Editorにあるのか，earthにあるのかは現在不明．CHAM社に問い合わせたところ以下のような回答があった．

ご指摘の件、はっきりとした原因は分かりませんが、現状で分かったことをご報告致
します。

slab1のy方向リージョンの分割数を10から9に変更するとこの隙間はなくなりまし
た。
ちなみにslab1のy方向リージョンの分割数を10から11とすると隙間が更に増えること
となりました。

toleranceの問題か確かめるために、
先生のｑ１をmm単位系で作り直してテストしてみましたが、この部分は分割数が10だ
とやはり
隙間ができ、９にすると解消されます。

ドメイン全体のサイズに対するリージョンのサイズ、
そのリージョンサイズの中で取り得る格子サイズの最小値がどこかで
決まってしまっているような気がします。

現状では，オブジェクトが境界を接するときには計算結果を注意して見るしか，この現象を防ぐ方法はないだろう．

ちなみに，オブジェクトどうしの隙間を発見した場合の対策は簡単．最小桁で「切り捨て」が出るようにしてやればよい．q1ファイルを手動で以下のように編集したところ，隙間は見事に消失した．

> OBJ, NAME, SLAB2
> OBJ, POSITION, 0.000000E+00, 0.000000E+00, 1.000000E-02
> OBJ, SIZE, 1.490000E-02, 1.501000E-04, 1.000000E-03

最終的に出来たq1ファイルは以下のようなものである．
q1ファイル


・計算結果

　計算結果は以下のようになった．


カーフの，切り口のごく近傍で圧力が±0.5atm近く変化しており，非圧縮近似がちょっと厳しいことを示唆している．流速は，切り口のところで最も速く400m/s程度．それ以外の，カーフ外側では300m/s程度とわかる．上，下ノズルの圧力はそれぞれ0.76kg/cm^2Gと1.64kg/cm^2Gであった．


C1，C2の分布から下ガスがカーフ近辺のガスの流れに与える影響がわかる．この条件では下ガスの勢いが足りず，ドロスフリー切断が期待できないことがわかる．

・まとめ