～ 目次 ～
1. やってみたこと
前回の記事で作成した洗面台の3Dモデルに水を満たし、パイプから水が流れ出ていくシミュレーションをやってみました。
2. シミュレーションの概要
洗面台に水を張り、水が重力で洗面台下部のパイプから流れ出ていくような解析を、自由表面解析（VOF法）でやっていきます。
そこで、incompressibleVoF の公式チュートリアル (waterChannel)1 をコピーして、メッシュや水の範囲などを変更することでシミュレーションモデルを構築しました。
OpenFOAM 12の場合は少し分かりにくいのですが、waterChannelのチュートリアルは、
$FOAM_TUTORIALS/incompressibleVoF/waterChannel/にあります。ここから、適当なディレクトリにコピーして変更を加えていきました。どこにどのような変更を加えれば良いか、初めは私も分からなかったのですが、Softflow(株)の「いきなりOpenFOAM」の連載2,3,4,5がすごく参考になりました。
3. シミュレーションモデルの構築
解析に使用した全ファイルはGitHubにあげておきました。使いたい方は、ダウンロードとかクローンするなりして使っていただければと思います。
今回はwaterChannelチュートリアルとの差分というか、変更の要点をピックアップして紹介したいと思います。
3.1. メッシュの作成
複雑形状のメッシュの頂点や面の定義を手動でやるのは骨が折れるので、Fusion 360（3DCAD）と XSim（プリプロセッサ）、snappyHexMesh（OpenFOAM内臓のメッシャー）を組み合わせてメッシュの作成を行いました。
3.1.1. Fusion 360 からサーフェイスを抽出
まずは、Fusion 360で3Dモデルからサーフェイスを抽出します（詳しいやり方はこちらの記事を参照）。抽出したサーフェイスをリネームして、side（排水タンク壁面と底面）, walls（液体が流れる洗面台の壁面）, atmosphere（洗面台の大気開放部分）, atmosphere2（排水タンクの大気開放面）として保存します。
walls以外の面は、Fusion 360のサーフェイス機能を使って作成しました（閉じた計算領域を作る必要があるため）。
3.1.2. XSim でメッシング設定
XSimのサイトにアクセスします。アクセスしたら、適当にプロジェクト名を入力して（ここでは「sink_sim」とかにしました）、「作成」ボタンを押すと新規プロジェクトが作成されます。
まずはスケールの変更を適用して、モデル1/1000倍にします。これで、OpenFOAMの世界でメートル(m)単位としてモデルを扱うことができます。
その次に、目標ベースメッシュ数をデフォルト値より2桁多くします。
水が接触する部分には3層のレイヤーメッシュを設定しました。
あとは物性、初期条件、境界条件などなど… 色々あるのですが、これらをすっ飛ばしてエクスポートします（メッシュの設定だけが欲しいため）。 他に必要な設定は後で行います。
今回使用しているOpenFOAMバージョンがOpenFOAM 12なので、フォーマットを「OpenFOAM 12」に設定してエクスポートします。
3.1.3. snappyHexMeshでメッシング
XSimからエクスポートしたフォルダを開きます。フォルダのルートディレクトリに移動したら、blockMesh → surfaceFeatures → snappyHexMesh の順番でコマンドを実行します。そうすると、constant/triSurface フォルダの中身に生成したメッシュが保存されるので、それをコピーして解析に使いました。
snappyHexMeshDictの中身は以下のように少し修正しました。XSimから出てきたものからほとんど変更していませんが、refinementSurfacesのlevelを(1, 1)に設定（このようにしないとParaViewでメッシュが見れなかったため）、locationInMeshの座標を変更しています（これはXSimでも設定できます）。また、regionsのtypeの設定については実際に存在する壁面は wallsとし、それ以外は patch に設定します。
ここで重要なのは locationInMesh の (x, y, z) 座標で、これが正しく設定できないとメッシングがうまくいきません…（汗汗）。
失敗したときの現象として、メッシュとメッシュの境界がうまく接続されない（例えば、atomosphereとwallsの接続部分のメッシュが歪んでいる）、Default Facesという謎の直方体メッシュが生成されるなどがあります。同じような現象になったら、計算領域が閉じていない、あるいは locationInMesh が間違えている可能性が高いので参考になさってください。あとは、失敗するときはメッシング時間がやけに長いです…(snappyHexMeshがうまくいくと150 [秒]くらいですが、 失敗すると400[秒]以上かかりました)
ケースディレクトリでparaFoamコマンドを打ち、paraViewを起動します。そうすると、メッシングうまくいっているかを視覚的に確認できます（OpenFOAMだと基本的に黒い画面とにらめっこなのでw）。
ちなみに先ほど話題にあげた locatioInMeshの座標ですが、paraViewの「Point Source」などを使うと確認できます。この「Point Source」の座標にlocationInMeshの座標を打ち込み、表示された点が計算領域の内部に入っていればOKです(‘ω’)ノ
system/snappyHexMeshDict の中身
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object snappyHexMeshDict;
}
castellatedMesh true;
snap true;
addLayers true;
geometry
{
atmosphere
{
type triSurfaceMesh;
file "atmosphere.stl";
regions
{
atmosphere
{
name atmosphere;
}
}
}
atmosphere2
{
type triSurfaceMesh;
file "atmosphere2.stl";
regions
{
atmosphere2
{
name atmosphere2;
}
}
}
side
{
type triSurfaceMesh;
file "side.stl";
regions
{
side
{
name side;
}
}
}
walls
{
type triSurfaceMesh;
file "walls.stl";
regions
{
walls
{
name walls;
}
}
}
}
castellatedMeshControls
{
features (
{
file "atmosphere.eMesh";
level 1;
}
{
file "atmosphere2.eMesh";
level 1;
}
{
file "side.eMesh";
level 1;
}
{
file "walls.eMesh";
level 1;
}
);
refinementSurfaces
{
atmosphere
{
level (1 1);
regions
{
atmosphere
{
level (1 1);
patchInfo
{
type patch;
}
}
}
}
atmosphere2
{
level (1 1);
regions
{
atmosphere2
{
level (1 1);
patchInfo
{
type patch;
}
}
}
}
side
{
level (1 1);
regions
{
side
{
level (1 1);
patchInfo
{
type wall;
}
}
}
}
walls
{
level (1 1);
regions
{
walls
{
level (1 1);
patchInfo
{
type wall;
}
}
}
}
}
refinementRegions {}
locationInMesh (0 0.03 0);
maxLocalCells 100000000;
maxGlobalCells 100000000;
minRefinementCells 1;
nCellsBetweenLevels 3;
resolveFeatureAngle 30;
allowFreeStandingZoneFaces false;
}
snapControls
{
nSolveIter 30;
nSmoothPatch 3;
tolerance 4.0;
nRelaxIter 5;
nFeatureSnapIter 10;
}
addLayersControls
{
layers
{
atmosphere
{
nSurfaceLayers 3;
}
atmosphere2
{
nSurfaceLayers 3;
}
side
{
nSurfaceLayers 3;
}
walls
{
nSurfaceLayers 3;
}
}
relativeSizes true;
expansionRatio 1.0;
finalLayerThickness 0.3;
minThickness 0.03;
nGrow 1;
featureAngle 60;
nRelaxIter 5;
nSmoothSurfaceNormals 1;
nSmoothNormals 3;
nSmoothThickness 10;
maxFaceThicknessRatio 0.5;
maxThicknessToMedialRatio 0.3;
minMedianAxisAngle 130;
nBufferCellsNoExtrude 0;
nLayerIter 50;
nRelaxedIter 20;
}
meshQualityControls
{
maxNonOrtho 65;
maxBoundarySkewness 20;
maxInternalSkewness 4;
maxConcave 80;
minFlatness 0.5;
minVol 1.00E-13;
minTetQuality -1e30;
minArea -1;
minTwist 0.05;
minDeterminant 0.001;
minFaceWeight 0.05;
minVolRatio 0.01;
minTriangleTwist -1;
nSmoothScale 4;
errorReduction 0.75;
relaxed
{
maxNonOrtho 180;
}
}
debug 0;
mmergeTolerance 1e-5;メッシングに成功すると defaultFaces は存在せず、下の画像のようなメッシュになります。
ちなみに失敗したときの画像はしたのような感じです↓
3.2. 材料物性
WaterChannelチュートリアルのデフォルト値を使いました。なので、粘性係数などもすべてデフォルト値です。
3.3. 境界条件
水の表面は大気開放にしていて、壁面は滑りなしの条件を与えています。計算領域全体に重力加速度が設定されているので、水が重力を受けて洗面台や下部のパイプに触れながら流れていく感じになります。
3.4. 初期条件（0ディレクトリの中身）
3.4.1. 流速U, 圧力p_rgh, nut, omegaなど
ほぼデフォルトのままなんですが、設定している境界面の名前がデフォルトと異なる場合は変更が必要です。
今回は walls, side, atomsphere, atmosphere2 という名前（これらはデフォルトと異なるものが混じっている）になっているので、ファイルを開いて名前や境界条件の種類などを修正していきます。境界条件はtype（静止壁、滑り壁、自然流入出など）やそのvalue（各境界条件の定義に必要な値）を適切に設定します。
3.4.2. 水の領域設定
今回は自由表面解析なので、計算領域のどの範囲が液体であるか指定する必要があります。
洗面台に水を満たすような感じにしたかったので、alpha.water.origファイルとsetFields中身は以下のようにしました。ちなみに「alpha.water = 1 だと全部水」、「alpha.water = 0.5 だと水と空気の境界」を表します。
SetFieldsDictで設定した直方体領域と交差する部分に液体部分が生成されます。設定する座標によって水面高さが決まるので、注意して設定していきます（私は桁を間違えて、水面高さがつるつるいっぱいな感じになりましたw）。
0/alpha.water.orig の中身
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object alpha.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [];
internalField uniform 0;
boundaryField
{
walls
{
type zeroGradient;
}
side
{
type zeroGradient;
}
atmosphere
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
atmosphere2
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
}
// ************************************************************************* //
system/setFieldsDict の中身
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defaultFieldValues
(
volScalarFieldValue alpha.water 0
);
regions
(
boxToCell
{
box (-0.200 -0.1375 -0.070) (0.200 0.1375 0.030);
fieldValues
(
volScalarFieldValue alpha.water 1
);
}
);
// ************************************************************************* //
3.5. 計算の設定
ほとんどwaterChannelのデフォルト設定を使っています。水は思いのほか高速で移動するので、時間の刻み幅deltaT を 0.001 [s]に、計算結果の書き出し間隔（writeInterval）を短くして 0.01 [s]にしました。あとは、洗面台の水は結構な量があったので、計算の終了時間 endTimeを 5 [s]に設定しています。
また、system/functionオブジェクトなどをコメントアウトしないとエラーになったので、コメントアウトしておきました（Geminiありがとう）。
system/controlDict の中身
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application foamRun;
solver incompressibleVoF;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5;
deltaT 0.001;
writeControl adjustableRunTime;
writeInterval 0.01;
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes;
maxCo 6;
maxAlphaCo 6;
maxDeltaT 1;
// ************************************************************************* //
4. シミュレーション結果
シミュレーション結果は下のGIF画像のようになりました。時間が経つにつれて、水がパイプを通って流れていくことが確認できます。
実際に同じような感じで洗面台に水を張って水を流してみたんですけど、もう少しゆっくり流れる印象でした（パイプの入り口に色々部品が付いているせいか、パイプの壁面が汚れとか製造方法の関係で管摩擦が発生するせいか、パイプ直径が実際はもう少し小さいせいか… 原因はいっぱいあるのでよくわからないです）
このシミュレーション結果から、ParaViewで水部分のメッシュを取り出して、Blenderでマテリアルとかライティングをしっかりやると下のようにリアルな画像が作れます！
参考：シミュレーション実行時間（42分）
-
使用したPCのスペック（ASUS ROG FLOW Z13, GZ301VF.316）
CPU：13th Gen Intel(R) Core(TM) i9-13900H基本速度: 2.60 GHz
ソケット: 1
コア: 14
論理プロセッサ数: 20
仮想化: 有効
L1 キャッシュ: 1.2 MB
L2 キャッシュ: 11.5 MB
L3 キャッシュ: 24.0 MB
※ 20コア（論理）並列計算をしました
![シミュレーション実行中のPC画面（CPU使用率100[%] !?）](pc_monitor_during_simulation.jpg)
5. 参考文献
-
GitHub(@OpenFOAM Foundation),「OpenFOAM-12/tutorials/incompressibleVoF/waterChannel」 https://github.com/OpenFOAM/OpenFOAM-12/tree/master/tutorials/incompressibleVoF/waterChannel ↩︎
-
SimFlow, いきなりOpenFOAM, 「逆さにした瓶からの液体の流出」, https://www.softflow.jp/tech-forum/ikinari-openfoam-49/ ↩︎
-
SimFlow, いきなりOpenFOAM, 「ボトルへの注水（その1）」, https://www.softflow.jp/tech-forum/ikinari-openfoam-42/ ↩︎
-
SimFlow, いきなりOpenFOAM, 「水路解析用のひな形をつくる」 , https://www.softflow.jp/tech-forum/ikinari-openfoam-35/ ↩︎
-
SimFlow, いきなりOpenFOAM,「液滴落下による波紋の解析」 , https://www.softflow.jp/tech-forum/ikinari-openfoam-47/ ↩︎