・第5回　実カーフ条件で仮想『水シミュレーション』-下準備

　前回までの作業で，水+空気，熱移動なしという条件なら液体がアシストガスで吹き飛ばされるケースの計算が可能となった．そこで，一気に実際の加工条件に行く前に，同様の条件で行われた実験とシミュレーションを比較して計算の妥当性を検証しておくのが良いだろう．

　下ガスノズルによるドロスフリー切断のメカニズムを解明する努力として，昨年(2002年)度は，CFDと『水シミュレーション』を行った．水シミュレーションとは，実カーフとおなじ形をアクリル樹脂で作り，流体金属の発生する位置に水を注射器で注入することにより，流体金属の挙動を模擬してこれを観察するものである．

これと同じものをPhoenicsで再現してみようと言うわけだ．いわば，シミュレーションのシミュレーション．

(1)形状ファイル作成
　前回までに作成したケースを元に，水シミュレーションのケースを作成する．まずはじめにやるべきは，実際の系を再現する形状ファイルの作成である．Phoenicsで形状ファイルを作成する方法は大きく4つある．他にもCADファイルをインポートするという方法があるが，高価な三次元CADが別に必要なのでここでは割愛．

  1. 標準で用意されている形状から選択
  2. Shape Makerを使う
  3. AC3D(CHAM版)を購入し，それを使う
  4. 自分でデータファイルを書く

今回作らなければならない形状は以下の通りで，それぞれ表のような手法で作成した．

形状データ 作成方法 備考
上ガスノズル Shape Maker
上ノズルのふた 標準形状ファイル hfcyl
上ガス流入境界 標準形状ファイル hfcyl
スラブ 標準形状ファイル cube14
水流入境界 標準形状ファイル default(cube3t)
流出境界 標準形状ファイル default(cube12t)
下ガスノズル 自分でデータファイルを書く データファイル作成スクリプト
下ノズルのふた 自分でデータファイルを書く 上に同じ
下ガス流入境界 自分でデータファイルを書く 上に同じ

カーフは長方形で，水は垂直の壁からにじみ出すという仮定を置いたが，実際の加工を模擬するときはより現実に近い方法に改める必要があると思われる．

　1～4の手法は必要に応じて以下のように使い分けている．まず1の手法．形状が円，矩形，円錐や，複雑だが自分で細部まで決定する必要のないもの(ダクト，自動車，椅子や机)などは，標準形状ファイルを使えば良い．

　2の手法は，それほど複雑ではないのに，ライブラリに見つからない形状ファイルを作りたいときに便利．Shape MakerはPhoenicsに同梱されて，たいていの場合はこれで用が足りる．たとえば上ノズルの様な形状には最も適した方法といえる．


　3の手法を活用したことはない．AC3Dは，3次元CGモデル作成の定番ソフトらしい．使ったことが無いのでわからないが，ほとんど想像できる限りのあらゆる形状が定義できるらしい．CHAM社は，AC3Dの供給元と交渉し，PhoenicsのVR Editorで直接使える形式(datファイル)が出力できるAC3Dを作らせて，Phoenicsのユーザーに有償配布している．価格は税別\10,000円(2003年現在)．本家のサイトから標準版が$49.95で買えること，　標準版AC3Dからdatファイルへのコンバートも可能なことは知っておいて損はない．

　最後の，自分で書く方法だが，datファイルの書式さえわかってしまえばこれが一番楽な場合が多い．とくに，形状が一連の数式で定義可能な場合は，適当なプログラムを作って，形状を微妙に変えたいときはパラメータを変えて実行させるという方法が一番効率がよい．私の場合は，多くのケースでこの方法を使っている．使用言語は，もっぱらPerl．いちいちコンパイルしなくても良いし，複雑な文字列操作もお手のもの．今回，下ガスノズルと下ガスの入り口を解析領域端部で閉じるためのプラグを出力するスクリプトはこんな感じ．出力結果は下の図のようになる．

(2)メッシュ作成
　続いて，これら形状ファイルを計算領域内に置いてゆく．VR Editorの操作性，レスポンスはこういうときにありがたい．ここで，一つ技を使っている．

　Phoenicsは構造格子を採用しているため，格子を等間隔に刻むと格子と形状の表面が(一般に)一致しない．これを避けるため，Phoenicsでは「リージョン」という概念が採用されている．これは，オブジェクトの端点を境界として計算領域を幾つかの「リージョン」に分割し，リージョンの境界にはかならず格子点が来るように格子分割が行われるというものだ．我々の問題にで，オブジェクトを素直に置いてゆくと下の図のようなリージョンが定義される．


　ところが，このリージョンというヤツは，ケースが複雑になるとかえって鬱陶しい場合がある．メッシュの数を増減したいときにはリージョン毎に独立して行う必要があるので，上の図の場合，x方向には10箇所のメッシュ数定義を行う必要がある．このケースでは，とくに下ノズルが形状定義領域いっぱいに存在しておらず，リージョン境界には実利的な意味もない．

　これを避けるために，オブジェクトがメッシュに影響を与えないように指定して，別の方法で明示的にリージョンを定義する方法がある．まず，オブジェクトのプロパティで"Object affects the grid"をNoに指定する．

これで，計算領域内のリージョンがぐっと少なくなる．スラブと空気の境界線のように，リージョンを生かしたい場合はそのままほおって置けばよい．

そして，必要な位置に明示的にリージョン境界を作るために，wirexyzオブジェクトを使う．このオブジェクトは透明な矩形で，物質としての実体は持たずにただその端点がリージョンを作る役割を持っている．これを使い，必要な位置にリージョンを確保したものがこちら．

次に，忘れてはいけないのが，上下の流入境界にリージョン境界が来るようにする操作．通常はオブジェクト境界に合わせてリージョンが作られるが，今回は手動で「オブジェクトをリージョン境界としない」ようにした．従って，流入境界の取り扱いはデリケートである．流入はz方向なので，x，yのメッシュがオブジェクトに合っている必要はないが，z方向には厳密にオブジェクト境界がメッシュの位置に来る必要がある．そのため，厚さ0.2mm，(x-y)には領域いっぱいのwirexyzオブジェクトを作り，それをz=0とz=0.02mの位置に置く．

あとは，必要なメッシュ密度を与えてやればできあがり．

(3)テストランその1：t=0からの過渡計算で水の落下を見る
　大規模な計算のばあい，計算が終わるまでに長いときでは数日かかる．最終結果を見て，実ははじめの仮定が間違っていたのに気がつくときほどむなしいものはない．したがって，計算も一気に目的のものを行うのではなく，簡単なものから慎重に進める必要がある．まずは，上ガス，下ガスノズルの流量を無視できるほどにして，t=0から0.1秒まで積分させてみた．液体の流出量は，今までのcut&tryから大きく変更しないように，0.1m/sとした．流量に換算すると0.15ml/sとなる．作ったq1ファイルは以下の通り．

解析領域 0.03×0.015×0.02[m]
メッシュ 80×80×80，不等間隔
Light fluid(気体) 窒素(@27℃)
Heavy fluid(液体) 水(@20℃)
境界条件 手前(対称境界)を除く全面流出境界
重力 -z，9.8m/s^2
スラブ(白) アルミニウム
液体流入量 0.15ml/s
上ガス流入量 1ml/s
下ガス流入量 1ml/s

q1ファイル


　
ご覧のように，t=0.070sで計算が発散している．原因としては，第一に空間のメッシュ密度が液体存在範囲で足りないということが計算結果から見て取れる．発散はちょうどwireオブジェクトのすぐ下，メッシュが粗になるところで始まっている．これは，液体に表面張力がないため，流れ出した液体が極細の糸になって流れ落ちるためである．したがってこれはメッシュ密度を上げれば解消する問題ではない．
　時間刻みについてはどうだろうか．自由落下する液滴が初速度0から1cm移動するのに要する時間は，空気抵抗を無視すれば44msである．したがって，アシストガスによる外力が重力と同程度，あるいは大きいレーザー加工の場合，決定的な現象はそれより短い時間で終了する．したがって，解析領域を1cm程度の大きさとすれば，44msより長い時間の積分は無意味とわかる．
　最後に，計算領域の大きさであるが，液体の動きとして本質的に重要なのは，上のケースでは中央の1/3程度にすぎない．比較的大きな解析領域を取ったのは，昨年度のアシストガスの計算と定量的な比較が可能になるためである．これも，メッシュ密度と計算速度のトレードオフを考えた場合，小さくした方が良いかもしれない．
　上の反省をふまえ，t=0でレーザーの真下に置かれた0.3mm×0.3mm×1mmの液滴の流れを追いかけることにした．時間刻みは0.5ms，積分は40msで打ち切る．周囲のメッシュ密度を下げ，カーフ近傍のメッシュ密度を上げた．

メッシュ 80×70×95，不等間隔
液体流入量 0.0

q1ファイル

　下の結果のように，重力によって落下する水が再現できた．

(4)テストランその2：アシストガスのみの定常計算
　液滴の落下は計算できたので，続いてアシストガスを導入する．ここで，今までは触れないでいたことを告白すると，SEM法による計算には，今回の目的に適合するかどうか微妙な，本質的な制限がある．それは，気体は非圧縮近似でしか解けないということである．SEM法では，液体と気体の密度の違いを利用して自由表面を表現しているので，密度が時間や場所の関数になっては都合が悪いというわけだ．一方，ノズルから吹き出すガスは1～2kg/cm^2Gなので，当然密度は倍以上変化する．また，流速も，ノズル出口で音速という縛りがあり，これがガスの運動量を支配する．これが非圧縮近似だとどうなるかを慎重に見極める必要がある．幸い，圧縮性流体で解いた同じ問題(自由表面はなし)があるので，これと比較することができる．

　まずは，テストラン1のケースを，「アシストガスのみの定常計算」に変形する．改修する点は

  1. Menu→Geometryで計算方法をTransientからSteadyに
  2. Menu→Numericsで緩和係数，Total number of iterationsの設定
  3. 水ブロックのオブジェクトを個体に変更．今回は「アルミニウム」とした．
  4. 上，下ガスの流入境界を設定．流速はとりあえずこのくらい．

これに加え，今回新たにわかった点についていくつか改修する．

  1. 液体ブロックを個体に変更しても，q1ファイルは自動更新されない．いったんVR Editorを閉じ，q1ファイルの初期条件記述エリアを手動で削除．
    Group 11.Initialise Var/Porosity Fields
    FIINIT(PRPS) = -1.000000E+00

    PATCH (LIQUID ,INIVAL,3,0,0,0,0,0,1,1)
    INIT(LIQUID ,SURN, 1.000000E+00, 1.000000E+00)
  2. SEMモデルにおいて，VFOL，SURNはWhole Fieldで解いてはいけない！！

    いままで，マルチノードの癖で解く(SOLVE)変数は全てWHOFに指定していた．これが，泥沼にはまった原因だった．SEMモデルにおける変数VFOLは，各パッチ毎に独立した値で，計算領域の中で保存される量ではない．WHOFを指定することにより，これを無理矢理保存するようにがんばった結果，INLET境界からガスが出てこないと言うトラブルが起こった．定常解析のときはSOLVE，STORE両方を切ったら上手くいったので，"WHOF=Y"が原因と気がつくまでに1週間を要した．

このようにして作ったq1ファイルを解いてみる．

オブジェクト テストラン1と同じ．ただし，液体ブロックは
同じ大きさのアルミブロックに置換
時間積分 steady
ガス流入条件 上ノズル：-5.0m/s
下ノズル：+5.0m/s

q1ファイル

　まあ，それなりに納得のいく結果を得た．乱流モデルも考慮していないし，ガス流量も適当なので，この段階ではまだ圧縮性流体の計算との比較はナンセンス．

(5)テストランその3：非定常解析のアタリをつける
　非定常解析の時間刻み幅は，流速の速いガスが存在するときは液体だけの場合に比べ非常に厳しくなると予測される．そこで，アタリをつけるために，テストラン1のケースにt=0からアシストガスを噴射する非定常計算をやってみた．時間刻みの目安は，離散積分のセオリーに従えば(最小の格子間距離)÷(ガスの流速)だが，上の例で考えると約1μsと，相当厳しい値となる．おそらく，観察したい現象は10msのオーダーで起こると思われるので，10,000枚ほどのタイムスライスを計算する必要がある．

オブジェクト テストラン1と同じ．ただし，アシストガスを噴射する．
時間積分 Transient Δtをいろいろ変える．計算枚数は1,000枚．
ガス流入条件 上ノズル：-5.0m/s
下ノズル：+5.0m/s

q1ファイル

　とにかく，上の様な流れ場で非定常計算が成立する時間刻みをしるため，Δtをいろいろと変えて計算を行ってみた．その結果，1000枚のタイムスライスを発散せずに完走することが出来たのはΔt=1μsとわかった．10ms後の流れ場を知るために10,000枚・・・・・しかも，ガスはまだ本気で吹いてない！

(6)テストランその4：定常解析(テストラン2)を初期条件とした非定常解析

オブジェクト テストラン3と同じ．ガス圧力，流速の初期条件は
テストラン2で得られた結果を使用
時間積分 Transient Δt=1μs．
ガス流入条件 上ノズル：-5.0m/s
下ノズル：+5.0m/s

q1ファイル

　ここで，テストラン2で作られたphyファイルを「rst」というファイル名にリネームする．これは，これから行う計算で作られるphyファイルが，初期条件のphyファイルを上書きしてしまわないためである．続いて，Initialzationパネルから，"Restart for all variables"をクリック．すると，初期条件が全て「READFI」に変わる(PRPSは例外)．
　しかし，定常計算と今回の計算は，液体の存在が異なるのでSURNの初期条件を引き継いではまずい．そこで，初期条件からSURNのところだけ，「1.00E-10」とする．最後に，リスタートファイルを「rst」としてできあがり．


　計算を実行してみたものの，面白い現象は何も見えなかった．ガスの流れははじめから最後までほぼ一定，液体も初期条件の位置から動くことは無かった．ほんの少し，表面がはがれた程度．やはり，積分時間が1msでは面白い現象を見るには短すぎる．

(6)テストランその5：時間間隔を長くする

計算条件 テストラン4と同じ．発散を防ぐため，SURNとVFOLの緩和係数を変える．
時間積分 Transient Δt=1μs．
ガス流入条件 上ノズル：-5.0m/s
下ノズル：+5.0m/s

q1ファイル

　上の計算では，水の挙動を見るために一週間もかかってしまう．そこで，Δtを短くすることを考えた．まず，Δｔを5μsに． 

 RSET(U,0.000000E+00,5.000000E-03,1000)
　このままでは計算がすぐ発散してしまうので，VFOLとSURNの緩和係数を変えた．
 RELAX(P1  ,LINRLX, 1.000000E-01)
 RELAX(U1  ,FALSDT, 1.000000E-03)
 RELAX(V1  ,FALSDT, 1.000000E-03)
 RELAX(W1  ,FALSDT, 1.000000E-03)
 RELAX(VFOL,LINRLX, 1.000000E-01)
 RELAX(SURN,LINRLX, 1.000000E-01)

　計算結果は以下のようになった．なお，この直後に計算は発散している．


　奇妙な計算結果だ．

　水の加速度は，水ブロック上にかかる圧力からほぼ概算できる．ブロック上部に掛かる圧力はt=250μsのとき約400Paと表示されている．

ブロックの断面積が0.15×0.3[mm^2]なので，ここにかかる力は約2×10^-5[N]となる．一方，ブロックの質量は0.15×0.3×1[mm^3]の水なので5×10^-8[kg]．加速度は約400m/s^2となる．はじめの3msでブロックが動くと思われる距離は1.8mm．悠々観察可能な大きさである．すなわち，上の計算結果は正しい物理現象を反映していないと言うことになる．

もういちど，基本に戻って計算をやり直さねば．

・まとめ