・第7回 SEMモデルに表面張力?副題:表面張力奮闘記

・プロローグ

 さて,前回の計算で,モデルに表面張力が考慮されていないため液滴の動きが非物理的になることがわかった.マニュアルによると,

RESTRICTIONS

The method cannot be used to produce steady-state solutions directly. However, steady-state conditions can be achieved asymptotically from a suitable transient simulation.

Surface tension effects at the interface are not included in the model. This is however a feasible extension to the model.

Numerical diffusion, despite treatments for its limitation, still imposes a restriction on minimum grid requirements.

Due to the explicit nature of the van Leer scheme, restrictions on the size of the time step increment can make long transient processes computationally expensive.

とあり,表面張力を独自に実装することは原理的に可能であると判断できる.実際,CHAMのホームページにも,表面張力を考慮した例題が幾つか出ている.しかし,これらの計算は標準パッケージに独自の実装を施したもので,これを実現するためにはPhoenicsソースコードへのアクセス権が必要だ.我々が購入したPhoenicsは,導入コスト削減のためアクセス権をつけていない.

一方,Phoenicsには,ソースコードにアクセスせずとも相当柔軟に流れの支配方程式を弄べる「in-form」という機能がある.マニュアルで調べ,「in-formでもでlきるのでは」と問い合わせたところ,

> いつもPHOENICSをご使用いただきありがとうございます。
>
> 1)informの可能性
> in-formでは界面形状から表面張力を計算するタイミングに処理
> を行えないので(タイムステップの最初か最後どちらでも良いのですが)
> 出来ません。groundファイルを使用するかgxsurfの中に組み込んで
> しまってください。
>
> 2)表面張力式
> 表面張力の計算は、色々な方法があるので教科書やローレンスリバモア
> のSOLA-VOFコードのマニュアルを参考にしていただければよいのですが
> 簡単な方法を紹介します。
> SEMをご使用であればVFOLが体積率なので曲率Kは
>
> K=-div((gradVFOL)/|grad(VFOL)|)
>
> で求まります。表面張力によるラプラス力は、表面張力係数σから
>
> f=σKdivVFOL
>
> になります。fを成分分解してu1,v1,w1のソースとしてに与えてください。

とのこと.groundファイル,gxsurfはソースコードのアクセス権が無ければ変更できない.絶体絶命か?しかし,ここでめげずに,マニュアルを調べて以下の方法では出来ないか問い合わせた.

> Example751のq1ファイルを見たところ,以下のような記述がありますが,
>
> ---------------------------------------------------------------------
>   Echo InForm settings for Group 13
>   INFORM13BEGIN
> PATCH(iREACRAT,VOLUME,1,NX,1,NY,1,1,1,1)
> (SOURCE of C1 at iREACRAT is RCONST*C1^EXPO*XG^(-2)*(1.0+YG)*(1-C1))
> (SOURCE of C2 at iREACRAT is RCONST*C1^EXPO*(1.0+YG)*(1-C1) WITH FI$
> XV)
>   INFORM13END
> ---------------------------------------------------------------------
> ここにU,V,Wのソースとしてf=σk dif VFOLを与えるという方法が考え
> られます(用意されている関数だけで実現するのは相当面倒ですが).
> 先ほどのメールによるとタイミング的にIn-formでの組み込み
> は無理とのご回答ですが,この方法をとった場合に発生するであろう
> 問題は何でしょうか.

すると,

> ご指摘のとおりSEMではVFOLはOLDの値が保持されているので、
> VFOLを使うのであればgroup13でも良いかもしれません。
> しかしinformではfの求解ができないので、fについて陽解法を使用
> することになりますので、timestepは相当細かくする必要がありそう
> です。さらにIZスラブごとにSURN->VFOLの代入を行いますので
> XY2次元の計算しかできません。

との解答.まあ,In-formに組み込むのだから,微分方程式に完全に組み込むのはどだい無理な話.しかし,SURNの値は毎ループごとに極端に変化するわけではないので,近似的にはいけるんじゃ無かろうかと見切り発車した.

そのあと,何回かのやりとりでIn-formに関する文法を教わり,In-formの中に表面張力を組み込んだ.


 ************************************************************
  Echo InForm settings for Group 13
  INFORM13BEGIN
STORE(SUR)
STORE(VFF)
STORE(GVX)
STORE(GVY)
STORE(GVZ)
STORE(GVA)
STORE(SUX)
STORE(SUY)
STORE(SUZ)
STORE(DVX)
STORE(DVY)
STORE(DVZ)
STORE(CVR)
PATCH(PA1,VOLUME,2,NX-1,2,NY-1,2,NZ-1,1,LSTEP)
(STORED of GVX at PA1 is (EAST(VFOL)-WEST(VFOL))/(DXG+DXG[-1]))
(STORED of GVY at PA1 is (NORTH(VFOL)-SOUTH(VFOL))/(DYG+DYG[&-1]))
(STORED of GVZ at PA1 is (HIGH(VFOL)-LOW(VFOL))/(DZG+DZG[&&-1]))
(STORED of GVA at PA1 is SQRT(GVX*GVX+GVY*GVY+GVZ*GVZ))
(STORED of SUX at PA1 is GVX/GVA with IF(GVA.GT.0.0))
(STORED of SUY at PA1 is GVY/GVA with IF(GVA.GT.0.0))
(STORED of SUZ at PA1 is GVZ/GVA with IF(GVA.GT.0.0))
(STORED of SUR at PA1 is SQRT(SUX*SUX+SUY*SUY+SUZ*SUZ))
(STORED of DVX at PA1 is (EAST(SUX)-WEST(SUX))/(DXG+DXG[-1]))
(STORED of DVY at PA1 is (NORTH(SUY)-SOUTH(SUY))/(DYG+DYG[&-1]))
(STORED of DVZ at PA1 is (HIGH(SUZ)-LOW(SUZ))/(DZG+DZG[&&-1]))
(STORED of CVR at PA1 is DVX+DVY+DVZ)
(SOURCE of U1  at PA1 is -CVR*GVX*VFOL*0.0740)
(SOURCE of V1  at PA1 is -CVR*GVY*VFOL*0.0740)
(SOURCE of W1  at PA1 is -CVR*GVZ*VFOL*0.0740)
(STORED of VFF is -CVR*GVA*VFOL*0.0740)
  INFORM13END
 ************************************************************

しかし,この形では計算が上手くいかない.

表面張力の計算で最もポピュラーなベンチマークは「振動する液滴」だろう.領域中央に球形の液滴を置き,上の式をIn-formで組み込むと,結果は下のようになる.

解析領域 0.01×0.01×0.01[m]
メッシュ 30×30×30,等間隔
時間刻み 50μs
Light fluid(気体) 空気(@20℃)
Heavy fluid(液体) 水(@20℃)
境界条件 全て対称境界(流出なし)
重力 なし
液滴のサイズ 6mmφの球体

q1ファイル

本来,計算結果は時間的に定常でなければならないが,計算を繰り返すうちに,変数VFOLがだんだん滲んできて,最後は液滴が解けて無くなってしまう.どうやら,In-formで力のsource項U1,V1,W1を与えていることが原因らしい.移流方程式の合間に,VFOLの滲みを解消する「再セット」操作をすればよいことは知られているが,これはそれこそソースコードにアクセスしないと無理な話.何とかならないかとCHAMに問い合わせたところ,

> <回答>
> 1)可能と思いますが流入流出がある場合、VFOLに関する境界条件を
>   変更する必要があります。
> 2)通常レベルセット法では何ステップかに1回体積保存を満足するように
>   距離関数の勾配がconstantになるようにニュートン法などで最初期化
>   をします。先生のpATCHの設定はおかしいです。DTはニュートン法の
>   仮想的なDTで在るべきです。
>   (最初期化しなければ距離関数も単にTVDで時間積分するだけなので
>    VFOLと同じ分布になります)
>   たぶんこれをやりたいのだと思いますがこれも時間ステップのタイミングで
>   VFOLを一気に修正する必要があるのでinformではできません。
> 3)これはユーザールーチンgxsurprp.htmで変更できればよいのですが
>   アクセス権が無いと修正できません。
> 4)よく考えたらこれも難しいです。前回
>     grad(H)=grad(DEN1)/(rho_liq-rho_gas)
>   にしてほしいといいましたが、これもおかしいです。
>     grad(H)=grad(DEN1)/(rho_liq-rho_gas) lolim<Vfol<uplim
>   が保証されないといけないのですがrho_liq,rho_gasが値を
>   持つため厳密には一致しません。
>
> 以上が貴質問と想定される件に関する回答です。
> アクセス件さえあればコーディング自身は難しいわけでは
> ありませんが、全てinformでは本来持っている機能の修正
> ですので、デフォルトのものをカットしたり、別の解法を
> 組み込んだりする必要などあり不可能だと思います。

いくつか,先方に事実誤認があったのでしつこく問い合わせたのだが,もううんざりらしく

> では本件のサポートは解決とさせていただきます。
> 今後ともよろしくお願いします。
>

と,勝手に「解決」にされてしまった.




・独自の実装

 CHAMのサポートは匙を投げたが,私の勘は,in-formの範囲で何とかなると告げている.まずは,in-form機能の勉強からスタートした.結構,様々な機能がある.特に,変数を計算するタイミングが制御できることがわかり,これが有効だと判明.

 ・・・それから2ヶ月,まあまあ納得の行く表面張力モデルを作り上げた.詳しい実装方法は企業秘密.

液滴の変形問題を解いたものを下に示す.液滴は予想通りまず菱形に変形し,そのあと頂点が押されて球形に戻る.どうだ!!

解析領域 0.01×0.01×0.01[m]
メッシュ 40×40×40,等間隔
時間刻み 50μs
Light fluid(気体) 空気(@20℃)
Heavy fluid(液体) 水(@20℃)
境界条件 全て対称境界(流出なし)
重力 なし
液滴のサイズ 4×4×4 [mm^3]



しかし,現状,問題はある.計算の安定度が不満足.本来,液滴は上下,左右対称に振動し,最終的には球形になって定常状態に落ち着くはずだが,■→◆→■くらいで,形状が余りはっきりしなくなる.何となく●になって,その後は不定形に振動する.

さらに,体積の保存が良くない.t=0から体積が徐々に増え始め,t=30msでは24%も増加する.この後体積は徐々に減り始め,t=0.1sで元に戻る.

これらの問題があることはあるが,in-formだけで表面張力を実装したのだから文句は言えないだろう.




・モデルの検証

 液滴の振動周波数は,液体の表面張力係数の関係で以下のように表される.参考リンク


従って,適当な初期条件を与えて液体を最低モードで振動させ,周波数を測れば表面張力モデルの妥当性が検証できる.

はじめは,初期条件を回転楕円体として最低次モードの振動を観測しようと思ったが,上述の不安定性により振動周波数を観測するのが無理だった.そこで,上の計算結果の振動周期から定性的にモデルの妥当性について言及するにとどめる.

上式より,4×4×4mmの液滴の最低モード振動周波数は16Hzとなる.上の計算結果の振動周波数は約20Hzだ.従って,定性的には計算結果は妥当であることがわかる.

まあ,この辺で勘弁してもらおうか.




・まとめ