重合格子とメニスカス壁
VOFで界面表現しようとすると、メニスカス壁が解析解より早くに壊れる
OpenFOAMの設定
モデルとしては円が徐々に下がっていくときにメニスカスがどの辺で破壊されるかを見たい。なので片方向連成でよい
背景メッシュのセル長：$0.2[mm]$
オーバーセットのセル長：$0.1[mm]$
snappyHexMeshのcastellatedMeshControlsは再分割を三次元に行うので、2次元計算をしたいのに円部分だけ三次元方向に複数メッシュを持ってしまう
そのため、extrudeMeshで奥行きのメッシュを1つにする
押し出す面とそれに相対する面を別個で指定するみたいなので、empty面はFrontとBackで分ける
oversetFront
{
type empty;
faces ((1 5 4 0));
}
oversetBack
{
type empty;
faces ((3 7 6 2));
}
重合格子側のsystemにextrudeMeshDictを作成する
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object extrudeMeshDict;
}
constructFrom patch;
sourceCase "./";
sourcePatches (oversetFront);
exposedPatchName oversetBack;
flipNormals true;
extrudeModel linearNormal;
nLayers 1;
expansionRatio 1.0;
linearNormalCoeffs
{
thickness 1e-3;
}
mergeFaces false;
mergeTol 0;
exposedPatchNameは、linearNormalでは必要性をあまり感じないが、選択しないとエラーを吐かれる
flipNormalsの順方向・逆方向の区別はワールド座標系基準。flipNormalsがtrueのときに逆方向（マイナス）なので注意
Allrunはこんな感じ
. $WM_PROJECT_DIR/bin/tools/RunFunctions
cd oversetmesh
blockMesh
snappyHexMesh -overwrite
extrudeMesh
cd ..
blockMesh
mergeMeshes . ./oversetmesh -overwrite
topoSet
restore0Dir
setFields
overInterDyMFoam
こんな感じになる
奥行以外の方向の再分割はそのままになっている
dynamicFvMeshはオーバーセットを使うのでdynamicOversetFvMesh
motionSolverは1つの壁を下に動かすだけなのでsolidBodyにする
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object dynamicMeshDict;
}
dynamicFvMesh dynamicOversetFvMesh;
motionSolverLibs ("libfvMotionSolvers.so");
solver solidBody;
solidBodyCoeffs
{
cellZone circleZone;
solidBodyMotionFunction linearMotion;
linearMotionCoeffs { velocity (0 0 -50e-3); }
}
linearMotionは並進
velocityは並進する速度ベクトル
単位は$m/s$
circleZoneはtopoSetDictでcellSetから定義する
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object topoSetDict;
}
actions
(
{
name background;
type cellSet;
action new;
source regionToCell;
insidePoints ((2.8e-3 0.5e-3 2.8e-3));
}
{
name circle;
type cellSet;
action new;
source cellToCell;
set background;
}
{
name circle;
type cellSet;
action invert;
}
{
name circleZone;
type cellZoneSet;
action new;
source setToCellZone;
set circle;
}
);
青がoriginal
赤はMesh2倍
四角いのは重合格子
2倍のメッシュ精度にすると曲率の正確性は増すが、ここで気になるのは破れている場所が異なること
originalはhole近傍で破れたが、2倍の方は背景側の水面が重合格子から離れると同時、スパッと切られたように破れた
重合格子内の曲率は背景の計算空間に広がっている水ありきなので、背景側の水面が重合格子の領域内にない場合、上に引っ張ろうとする力を収束計算できない
バンジージャンプの紐を切られたみたいな状態
改善するには、重合格子を上に伸ばせばいいのかと思いきや実は上だけでは足りない
解析解超えちゃってるし、下図のように界面もおかしい
つまりは、重合格子の横幅が狭すぎて背景側への収束計算の影響力が小さい
→先んじて背景側のメニスカス破壊が起きて自由水面に戻ろうとする
→それに引っ張られて重合格子側の水面曲率が小さくなる
→メニスカスが内側に閉まっていく力の抑制になり、破れなくなる
こと界面においては、重合格子の包含する領域が計算の正確性にかなり影響することが分かる
横にも広げてみた
（下図は界面の曲率計算にLaplacianSmoothingとDensityWeightedを使っている）
さらに横を倍にした
メッシュ解像度は同じだが、グレーディング使って円近傍のセル長を先のケースと同じにしている
結果としてはほぼ変わらん
これ以上伸ばしてもあまり意味はないようだ
グレーディングで円近傍のメッシュを細かくした
比較してかなり良い挙動になった（破壊後の水面は怪しいが）
解析解よりも沈んでいるように見えるが、境界条件の問題と思われる
破壊が起こっている場所はほぼ解析解と同じである
円の境界条件はmovingVelocityWallで、要は速度を持ったnoSlip境界条件
このnoSlipが、ミクロの世界ではあまり良くないとされている（今回のケースではどうだか不明だが）
何はともあれ、メニスカス壁の表現には少なくともこれくらいの重合格子は必要であることがわかった。
オーバーセット内の横のメッシュ数を2倍にしたら、伸びてしまった
原因としては、水面の曲率がなかなか変化しないためだと思われる
円の境界条件、計算空間の境界条件（左右）の問題も考えられるが、どうやら密度重みづけが悪いらしい
消したら治った
結局、横のメッシュを細かくすると若干破壊のタイミングが早くなった
メッシュは正方形の方が安定するのか、円表面の滑りか、はたまたそのどちらもか…
とはいえ、ある程度の精度が2Dで実現できたので、今度は3Dで解析を行うことにする
