研究メモ ~SpuriousCurrent編~
SpuriousCurrentとは、VOF法で表面張力のエラー（誤差）によって界面近傍にありえない速度が発生する現象のこと
VOF法はセルごとの体積率によって液相を定義するため、表面張力は体積率をもとに計算される
$$ F_{\sigma}=\sigma k \mathbf{n}\delta_{s} $$
$F_{\sigma}$：表面張力
$\sigma$：表面張力係数
$k$：曲率
$\mathbf{n}$：界面単位法線ベクトル
$\delta_{s}$：デルタ関数
ソースコードを見てみる
interfaceProperties.C（106~147）
void Foam::interfaceProperties::calculateK()
{
const fvMesh& mesh = alpha1_.mesh();
const surfaceVectorField& Sf = mesh.Sf();
// Cell gradient of alpha
const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
// Interpolated face-gradient of alpha
surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
//gradAlphaf -=
// (mesh.Sf()/mesh.magSf())
// *(fvc::snGrad(alpha1_) - (mesh.Sf() & radAlphaf)/mesh.magSf());
// Face unit interface normal
surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
// surfaceVectorField nHatfv
// (
// (gradAlphaf + deltaN_*vector(0, 0, 1)
// *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
// );
correctContactAngle(nHatfv.boundaryFieldRef(),gradAlphaf.boundaryField());
// Face unit interface normal flux
nHatf_ = nHatfv & Sf;
// Simple expression for curvature
K_ = -fvc::div(nHatf_);
// Complex expression for curvature.
// Correction is formally zero but numerically non-zero.
/*
volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
forAll(nHat.boundaryField(), patchi)
{
nHat.boundaryFieldRef()[patchi] = nHatfv.boundaryField()[patchi];
}
K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
*/
}
これは曲率$k$の計算
const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
→セル中心の体積分率（スカラ場）の勾配
つまりセルの中心から伸びる、水面に対する法線ベクトル
surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
→体積分率の勾配をセル中心から界面中心に変換する
surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
→界面フェイス単位の界面法線ベクトルを計算
ここで、
deltaN_
(
"deltaN",
1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
),
$\delta_{N}$は$|(\nabla\alpha)_{f}|$が限りなくゼロに近いとき、数値発散を防ぐ安定化項
（体積分率0.00…001のセルの界面法線ベクトルとか）
$\sum\frac{V_{i}}{N}$は全セル体積の平均
立方体だとすると、セル長$x$の3乗が体積なので、1/3乗してセル長の平均となる
即ち、
平均セル長が$1\mathrm{m}$のときに$10^{-8}$だけ分母に足されるのを基準として、
全体的にセルが小さいほど$\delta_{N}$は大きくなり、安定性を上げようとする
全体的にセルが大きいとその逆
セルがとても小さいときに0.9999…みたいな体積分率ができやすくなるから？
alpha1は第1相の体積分率らしい（第1相って気相液相どっち？？）
ちなみにdeltaN_のコードはOpenFOAM-5.0とかのinterfaceProperties.Cにあり、ESI版(OpenFOAMv2012)だと以下のようになっている
deltaN_
(
"deltaN",
1e-8/cbrt(average(alpha1.mesh().V()))
),
1/3が無くなっていることから察せるように、cbrtは立方根（cube root）のことである
cbrtの処理はdimensionedScalar.Cにある
dimensionedScalar cbrt(const dimensionedScalar& ds)
{
return dimensionedScalar
(
"cbrt(" + ds.name() + ')',
pow(ds.dimensions(), dimensionedScalar("(1|3)", dimless, 1.0/3.0)),
::cbrt(ds.value())
);
}
余談おわり
こんな感じで作った
メッシュは以下の通り
blocks
(
hex (0 1 2 3 4 5 6 7) (100 1 100)
simpleGrading (((0.35 0.2 0.3) (0.3 0.6 1) (0.35 0.2 3.3)) 1 ((0.35 0.2 0.3) (0.3 0.6 1) (0.35 0.2 3.3)))
);
最大時間刻みは1e-8で、クーラン数0.2
なんか思ってたんと違う
文献では対称的な流れ場が形成されるのに対し、こちらでは上部にだけ渦ができる
しかし偽流のメッシュ依存性は確認できた
時間ごとの最大流速
とんでもない変わりよう
あまりにも差がありすぎて極端な例なんじゃないの？って思ったので
研究で使っているアメンボモデルと同じスケールで試してみた
絶賛解析中の2jointモデルをcheckMeshすると、
Checking geometry...
Overall domain bounding box (-0.041 -0.041 -0.03) (0.041 0.041 0.03)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (2.3132359e-17 -6.923056e-18 2.8422549e-17) OK.
Max cell openness = 3.2930907e-16 OK.
Max aspect ratio = 6.4966365 OK.
Minimum face area = 1.8076838e-10. Maximum face area = 2.0769677e-06. Face area magnitudes OK.
Min volume = 3.7027244e-15. Max volume = 2.9461787e-09. Total volume = 0.00040445564. Cell volumes OK.
Mesh non-orthogonality Max: 67.722716 average: 5.3396825
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.5316666 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
Minimum face area（$= 1.8076838\times10^{-10}$）からセル長を求めると、
$\Delta x_f= ( 1.8076838\times10^{-10})^{\frac{1}{2}}\approx 1.3445\times 10^{-5}[\mathrm{m}]$
Min volume （$= 3.7027244\times10^{-15}$）からセル長を求めると、
$\Delta x_v= (3.7027244\times10^{-15})^{\frac{1}{3}}\approx 1.54706\times 10^{-5}[\mathrm{m}]$
大体$13[\mathrm{\mu m}]\sim15[\mathrm{\mu m}]$くらい
2次元液滴ケースの一番ちっちゃいセル長はsimpleGradingで中心部のセルなので、
ここで、$L$は計算空間の1辺長、$C$は計算空間1辺の分割数で、あんまり計算コストを上げたくないため、$C=100$とすると、
というわけで、simpleGradingは変えずに$3\times3[\mathrm{mm}]$の計算空間で分割数を$100\times100$にすればよい
アメンボのケースに合わせて、時間刻み1e-6、最大時間刻み1.96e-6、クーラン数0.6
2jointモデルのケースでは、最大時間刻みを次式をもとに設定している
いろんな論文で採用されてるやつ
$\rho$は液体の密度
本来は2相の密度の和だが、空気は十分軽いとして密度0で考えている
この式の導出についてはここで説明した（正しいかは知らん）
できた
解析してみる
やはりオーダーを変えてもメッシュ依存性は顕著なようだ
肢を太くしたら計算が発散したの原因は、おそらくこのメッシュ依存のSpuriousCurrentでほぼ間違いないだろう
肢周りの表面張力が加速→速度依存のダンパが爆発
S-CLSVOF solver
S-CLSVOFのソルバライブラリが配布されていたので使おうとしたが、ESI版の対応バージョンがv1812とv2212のみ
弊環境のv2012でv1812のソルバをwmakeするとエラーを吐かれるため
LibrunとLibcleanを作った
Allwmake、Allwcleanの代わりに実行する
本質的にはビルド順を変えて依存関係を解消した
Librun
#!/bin/sh
cd interfacePropertiesPsi
wmake libso
cd ..
cd twoPhasePropertiesPsi
wmake libso
cd ..
cd immiscibleIncompressibleTwoPhaseMixturePsi
wmake libso
cd ..
wmake 2> wmake.log
Libclean
#!/bin/sh
cd immiscibleIncompressibleTwoPhaseMixturePsi
wclean
rm -r lnInclude
cd ..
cd interfacePropertiesPsi
wclean
rm -r lnInclude
cd ..
cd twoPhasePropertiesPsi
wclean
rm -r lnInclude
cd ..
wclean
rm $FOAM_USER_LIBBIN/libimmiscibleIncompressibleTwoPhaseMixturePsi.so
rm $FOAM_USER_LIBBIN/libinterfacePropertiesPsi.so
rm $FOAM_USER_LIBBIN/libtwoPhasePropertiesPsi.so
しかしwmakeは通らない
ルートディレクトリのwmakeでエラーが起きるみたい
alphaEqn.H:136:27: error: no matching function for call to ‘correct(Foam::volScalarField&, Foam::surfaceScalarField&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, int, int)’
136 | MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
| ~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
OpenFOAMv1812とOpenFOAMv2012のalphaEqn.Hを比較すると、
OpenFOAMv1812
MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
OpenFOAMv2012
MULES::correct
(
geometricOneField(),
alpha1,
alphaPhi10,
talphaPhi1Corr0.ref(),
oneField(),
zeroField()
);
MULESの仕様が変わったみたい
alphaEqn.Hの中に、MULES::correctが2か所、MULES::explicitSolveが1か所あるのでOpenFOAMv2012の記法に書き換えてみると、wmakeが成功した
ディクショナリファイルの変更
0.origの中にpsiを追加する
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1812 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object psi;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
wall
{
type zeroGradient;
}
front
{
type empty;
}
back
{
type empty;
}
}
// ************************************************************************* //
system/fvSolution内のsolversの中に以下を追加
psi
{
$alpha.water;
nAlphaCorr 3;
nAlphaSubCycles 2;
cAlpha 1;
icAlpha 0;
scAlpha 1;
deltaX 1.5e-05;
gammaCoeff 0.75;
epsilonCoeff 1.5;
deltaTauCoeff 0.1;
}
deltaXは界面近傍のセル長で、先で計算したように、$0.015[\mathrm{mm}]$に設定している
gammaCoeff, epsilonCoeff, deltaTauCoeffは良くわからんのでtutorialの数値のまま
それ以外はalpha.waterと同じ
controlDictのapplicationを変える
application sclsVOFFoam;
結果
見た目ではいい感じになってるけど、グラフで見ると実は偽流は収まっていない
S-CLSVOF solverは、CSFモデルで計算された表面張力の向きのみを補正するソルバである
つまりこの結果で正しく動作していることが分かる
interFoamでは、表面張力誤差による偽流に押されて液滴が下に動いていくが、sclsVOFFoamではうまく中心に留まっている
表面張力の向きという観点ではinterFoamよりも物理的に正しく、S-CLSVOFの目的は
果たされているのだ
このソルバの欠点としては、壁境界条件が使えないことだろうか
接触角の設定ができないらしい（実装されていないだけ）
ていうかoverInterDyMFoamが使えなければ意味がないので、この解析は要らんかったかもしれん
Laplacian-smoothed
例のごとく、v2012用にソースコードを書き変える
interfaceProperties.C
convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)
↓
degToRad() * acap.theta(U_.boundaryField()[patchi], nHatp)
ラジアン変換の関数が違うみたい
この関数の定義元をincludeしておく
#include "unitConversion.H"
ちなみにこの関数は引数無しだと、ラジアン変換するときにかけるやつ
//- Multiplication factor for degrees to radians conversion
inline constexpr scalar degToRad() noexcept
{
return (M_PI/180.0);
}
こいつは要らないので消す
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
const Foam::scalar Foam::interfaceProperties::convertToRad =
Foam::constant::mathematical::pi/180.0;
alphaContactAngleFvPatchScalarField.Hは後のバージョンで静的、動的接触角に分かれたので、
全て`constantAlphaContactAngleFvPatchScalarField.Hに書き換える
cAlpha_
(
readScalar
(
alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
)
),
↓
cAlpha_
(
alpha1.mesh().solverDict(alpha1.name()).get<scalar>("cAlpha")
),
rhoとかphasesは型が分からないのでlookupのままでいいや
interfaceProperties.H
interfaceProperties(const interfaceProperties&);
void operator=(const interfaceProperties&);
↓
//- No copy construct
interfaceProperties(const interfaceProperties&) = delete;
//- No copy assignment
void operator=(const interfaceProperties&) = delete;
これはもう使わないので消す
//- Conversion factor for degrees into radians
static const scalar convertToRad;
Make/options
-I$(LIB_SRC)/transportModels/twoPhaseProperties/alphaContactAngle/alphaContactAngle \
↓
-I$(LIB_SRC)/transportModels/twoPhaseProperties/alphaContactAngle/constantAlphaContactAngle \
こいつを追加
-I$(LIB_SRC)/OpenFOAM/global/constants
/curvatureModel/vofsmooth/vofsmooth.C
財団版独自の型名が使用されているので、ESI用に置き換える
How to implement reconCentral in Openfoam1812
const unallocLabelList& pFaceCells = mesh.boundary()[patchi].faceCells();
↓
const labelUList& pFaceCells = mesh.boundary()[patchi].faceCells();
下図の"laplacian 2"はLaplacian-smoothedを1ステップで2回かけたもの。“laplacian 2 DW"はそれに加えて密度重みづけを付与したもの
結果的には、“laplacian 2 DW"が最も効果的なように見えるが、界面の曲線が平滑化されすぎるので、使わない方が良い
そのため、最も有効だと考えられるのは"laplacian 2"で、かつoverInterDyMFoamでも使用できる
この結果からは表面張力の正確性にどれだけ影響するのかわからない。あくまで表面張力から生まれる気相の偽流が軽減したことを示している
