・第1回 まずは,液体を落としてみよう

 未知の問題を解こうとするとき,先ず見るべきはEXAMPLESである.ところが,今回は,三次元,自由表面の良い例題が見つからなかった.「Library→Option Library→Advanced multi-phase→Equation Model(SEM)と辿ると例題がある」とマニュアルにはあるが,二次元の問題ばっかり.本質的に時間の掛かる三次元の問題は,例題には向かないのだろうか.

 そこで考えたのが,以下のような系である.領域は0.3m×0.3m×1mで,一番上に水の固まりを置く.t=0でこれが自由落下を始めたときに,その動きを追いかけようというものだ.東西南北の境界は対称境界として,下に流出境界を置く.


Falling water blocの計算

なぜ,水が流入する問題から始めなかったかというと,どうやっても上手くいかなかったから.あまりにくだらない袋小路は,このページでは切ることにしよう.

 以下の説明は,Phoenicsを使ったことが無い人には具体的なイメージがわかないかもしれないが,その場合は読み飛ばしてほしい.


ケースの作成手順

  1.  VR-Editorから,上のような大きさの領域を作成,ブロックと流出境界を一つずつ定義する.ブロックの大きさは0.1m×0.1m×0.1mで,上端にぴたりとつけた.
  2. Geometryの設定は以下の通り.このページで設定したのは,SteadyをTransientに変えたこと,そしてNumber of cellsをx,y,z各20に定義したことだ.
  3. 続いて,Time step settingsの設定に入る.これは,「計算は何秒間にわたって行い,それをどの程度のステップで刻むか」を設定するものである.入力するのは,当面は上の4つのボックスのみである.この例では,0から1sまでを500ステップに割っている.1ステップ当たり,2ms進む計算である.
  4. Modelの設定.ここで,当面は乱流は考えないことにしよう.設定するのは,Free-surface modelsをSCALAR EQUATIONにすることと,Solution for velocities and pressureをonにすることである.
  5. ここで,試行錯誤から得たノウハウに言及しておく必要がある.PhoenicsのSEMモデルをParallelで計算させるときは,Solution controlをいじらないと計算が始まらない.P1,U1,V1,W1,VFOL,SURN(要するに,解く変数全て)のWHOFを"Y"にしなくてはいけない.この関門は,超音速流れを解くときにCHAM Japanから教わったので比較的楽に通過出来た.

    03/09/04追記 CHAM-Japanサポート担当より,驚愕の事実を告げられる.SEMモデルはParallel版Phoenicsに対応していないのだそうだ.この前のセミナーではそんな話は聞かなかったのだが(怒).したがって,今までの計算は,偶然動いたのか,またはSEMはParallelで一応計算できるが,CHAMとしては正式にサポートしていないかのどちらかである.いずれにしても,Parallelで計算することは論外.このシリーズの計算は全てシングルノードでやり直した.WHOFはデフォルトでは"N"であるが,シングルノードの場合には実害は無いはず.
  6. Propertiesはいじらない.デフォルトの設定は,「軽い流体は空気,思い流体は水,t=0で系は軽い流体で満たされている」条件である.
  7. 続いて,Sourcesをいじらなければいけない.今回の計算は,重力で落ちる水の固まりだから,外力として重力(浮力)項を与える必要がある.設定は簡単,簡単.浮力のモデルにはCONSTANT,Density Difference,Boussinesq,Linear in 2 scalesがあるが,一番単純なCONSTANTを選択する.重力の方向は-zとする.後の設定はデフォルト.
  8. Numericsの設定.Phoenicsの実時間解析の考え方は次のようなものである.「ある時刻において,収束が得られるまで繰り返し計算を行い,収束が得られるか,設定されたループ回数に達したら時刻をΔt進める」.従って,定常解析と同様,Numericsの設定は,収束を得るために非常にクリティカルなものである(はずだ).しかし,一つ一つのタイムスライスで1,000ループも回していたのでは,計算が永久に終わらない.このケースではTotal number of iterationsを100に設定し,Relaxationはデフォルト値で一切いじらないようにした.これでも,各タイムスライスで10^0程度の収束は得られる.
  9. 最後に,Outputの設定.これをしないと途中経過が出ない.Field dumpingをクリックして,「どれくらいの割合で途中経過を出すか」指定する.Frequencyを20とすると,500ループが25枚のスナップショットになる.まあ,これくらいでよいだろう.
  10. さて,計算条件の設定のあとは境界条件だ.設定するのはブロックの物性と流出境界.まず,流出境界から行こう.ここでは,まずは何も考えずに,External valuesを全てIn-cellに設定しておく.External pressureはもちろんデフォルト値.このあたりの設定は,より現実に近い計算ではいじる必要が出てくるだろう.
  11. 水ブロックの設定.これが少々ややこしい.まず,オブジェクトのatrributeでTypeをUSER_DEFINEDに設定する.

    ついで,Attributesは以下のようにする.

    図では,PATCH Numberは1/1になっているが,これは一つ新しいものを定義したから.やり方は以下の通り.
    1. Newのところに適当な名前(FLUID)を入力.いったんOKをクリック.もう一度,Object Atrributeでこの画面に戻ると,上のようにPATCH number1/1がセットされる.
    2. 続いて,Name FLUIDの隣,Typeのところに"INIVAL"と入力.これは,このブロックに設定されるのは初期条件という属性ですよと宣言することに相当する.ITFとITLはどちらも1にセット.これで,t=0以外の時刻ではここに明示的に水の存在が定義されることはない.
    3. 続いて,VariableのSURNを1.0に,Coefficientを1.0にする.SURNは水であることを示す変数で,0~1の値を取る.Coefficientの1.0は,おまじないと思っておいて下さい.

 これで,準備完了.出来たものがこちらにございます.

 さて,二つの計算結果をお見せしよう.上の図が,いま紹介した方法で作ったもの,下の図が,time stepだけを10msとしたものだ.変数SURNは,上で述べたように水の存在を表す変数である.
 Time stepが粗い方は,途中から水が霧のように雲散霧消してしまう.もちろん,上の図の方がより現実の世界に近い訳だが,下のケースが上手くいかなかった理由は,時間方向の刻みの荒さだ.Transientな計算が上手くいくためには,Δt離れた二つの解析領域で,流体が大きく動いてはいけないと言うことである.今回の計算では,一番下にいるときの速さがおよそ4.4m/s,固まりの大きさが10cmだから,固まり一つ分動くのに最小20msしかかからない.また,SURNのAverage valueを見てゆくと,上の結果でも,水の質量が保存されているとは言えないことがわかる.これは,空間方向の刻みの荒さが原因と思われる.
 上の図の結果が正しいという保証もないのだが.水が落下しながら空気抵抗で広がってゆく様子は,現実にもこのような現象が見える.系は横方向に対称境界としているので,実際には20cm×20cm×10cmのブロックを落としていることに注意.空気抵抗によって中央が押し広げられ,ドーナツの様になって落下する様子が再現されている.更に,中央部が下端に到達する時刻も,重力加速度から単純に計算した値が0.44sだから,妥当な線だろう.
Result1

Result2
 一方,Numericsの繰り返しをどこで打ち切るかも,計算条件最適化に重要なファクターと考えられる.ソルバーの実行画面を見ていると,各時刻およそ30ループくらいで収束に到達しているように見えるが,敢えて繰り返しを2でカットしてみよう.
Result3
驚いたことに,計算結果はほとんど変わらない.ソルバーの画面を見ていても,最初の2ループで劇的にエラーが減少するので,100ループはこのばあい多過ぎということがわかる.これは,Transientの「収束」と,Steadyの「収束」の意味が違うということに起因する.Transientの場合は,たとえ液体が空中に浮いていても,その瞬間に(前time sliceからの)保存則が満足されていれば良いので,拘束条件は定常状態に比べて遙かにゆるい.一方,定常解を得るための繰り返し計算は,各格子点においてdφ/dt=0(φはあらゆる変数)が得られるまで計算を続ける必要があり,緩和係数を適当に設定しないと「振動」を引き起こすことになる.これが,Transientのケースで繰り返しが本質的に重要でない理由である.当たりをつけるために,繰り返しを2で止める方法は有効だろう.


・まとめ

03/09/04追記 Result1の計算結果をSingle node版で計算し,Parallel版と比較してみた.その結果,下図の用に両者の結果は全く一致.どうやら,SEMモデルでパラレルが使えないとは必ずしも言い切れないようだ.
Parallel Phoenics(Result1)Single Phoenics