はじめに
※ 今回は、既存の情報を自分のためにまとめただけです。先達はあらまほしきかな。深く感謝申し上げます。
早速ですが、昨年末にOpenFOAMで水槽の餌やフンなどの残渣を漉しとる装置を開発してほしいという案件がありました。
想定している水中浮遊残渣自体は、若干水よりも重いので沈むのと、水を含んで柔らかいので、単純に水中ポンプなりマグネットポンプで吸い出してやれば良さそうです。
ただし、そのまま捨てると水がもったいないので、再処理して（オゾン殺菌とかだ思いますが）水槽に戻したいということだそうです。あと、残渣を漉しとって乾燥させて肥料に使えないかということで、ちょっとシミュレーションをもとに設計してみようかということになりました（同じ学内の先生からの相談で、お金の話には未だなってません。ほんと自由にトライできるだけの研究費がないなと思います。校費があっても、ほとんど学生の材料費で消えます。）。
OpenFOAMであれば、解析できそうなので、対応するソルバーをいろいろ調べたところ、Penguintis様のサイトに詳しくまとまっていました（バージョンが古いけど、本家よりも分かりやすい）。
http://penguinitis.g1.xrea.com/study/OpenFOAM/particle_solver/particle_solver.html
- 熱計算不要
- 粒子同士の反応は不要
- 2相流
という仮定を行って、さらに調べたことろ、ESI版v1812にMPPICInterFoamというちょうど良さそうなソルバーを発見し、これをを使ってみることにしました。ちなみにOpenFOAM6には、MPPICFoamとMPPICDyMFoamだけです。
このほかDPMFoamというソルバーもありますが、DPM系は粒子密度（正確にはパーセル密度か？）が多い場合には精度が良いそうですが、計算時間はMPPIC系よりも大幅にふえるとのことです。
MPPICInterFoamとInterFoamの違い
ソルバーのソースコード(いずれもv1812のもの)のdiffをとってみると以下のようになりました。
MRF(fvOption)は使えるようですが、dynamicMesh機能は使えなさそうです。
Interfoamに、Multi-Phase Particle In Cell (MP-PIC) Lagrangian の機能を追加したものらしい。
MP-PICは計算効率が良く、いろいろ汎用性はありそうです。
（なんでFoundation版にはないのかなと思いました。）
$ diff MPPICInterFoam.C ../interFoam/interFoam.C
5c5
< \\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
---
> \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
25c25,28
< MPPICInterFoam
---
> interFoam
>
> Group
> grpMultiphaseSolvers
29,35c32,34
< (volume of fluid) phase-fraction based interface capturing approach.
< The momentum and other fluid properties are of the "mixture" and a single
< momentum equation is solved.
<
< It includes MRF and an MPPIC cloud.
<
< Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
---
> (volume of fluid) phase-fraction based interface capturing approach,
> with optional mesh motion and mesh topology changes including adaptive
> re-meshing.
39a39
> #include "dynamicFvMesh.H"
45d44
<
47c46
< #include "PhaseCompressibleTurbulenceModel.H"
---
> #include "turbulentTransportModel.H"
53,54d51
< #include "basicKinematicMPPICCloud.H"
<
61,63c58,61
< "Solver for two incompressible, isothermal immiscible fluids using"
< " VOF phase-fraction based interface capturing.\n"
< "Includes MRF and an MPPIC cloud."
---
> "Solver for two incompressible, isothermal immiscible fluids"
> " using VOF phase-fraction based interface capturing.\n"
> "With optional mesh motion and mesh topology changes including"
> " adaptive re-meshing."
71,73c69
< #include "createMesh.H"
< #include "createControl.H"
< #include "createTimeControls.H"
---
> #include "createDynamicFvMesh.H"
74a71
> #include "createDyMControls.H"
77,78c74,75
< #include "createFvOptions.H"
< #include "correctPhi.H"
---
> #include "initCorrectPhi.H"
> #include "createUfIfPresent.H"
84d80
< #include "readTimeControls.H"
90d85
<
95c90
< #include "readTimeControls.H"
---
> #include "readDyMControls.H"
112,149d106
< Info<< "Evolving " << kinematicCloud.name() << endl;
<
< kinematicCloud.evolve();
<
< // Update continuous phase volume fraction field
< alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
< alphac.correctBoundaryConditions();
<
< Info<< "Continous phase-1 volume fraction = "
< << alphac.weightedAverage(mesh.Vsc()).value()
< << " Min(alphac) = " << min(alphac).value()
< << " Max(alphac) = " << max(alphac).value()
< << endl;
<
< alphacf = fvc::interpolate(alphac);
< alphaRhoPhic = alphacf*rhoPhi;
< alphaPhic = alphacf*phi;
< alphacRho = alphac*rho;
<
< fvVectorMatrix cloudSU(kinematicCloud.SU(U));
< volVectorField cloudVolSUSu
< (
< IOobject
< (
< "cloudVolSUSu",
< runTime.timeName(),
< mesh
< ),
< mesh,
< dimensionedVector(cloudSU.dimensions()/dimVolume, Zero),
< zeroGradientFvPatchVectorField::typeName
< );
<
< cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
< cloudVolSUSu.correctBoundaryConditions();
<
< cloudSU.source() = vector::zero;
<
152a110,148
> if (pimple.firstIter() || moveMeshOuterCorrectors)
> {
> mesh.update();
>
> if (mesh.changing())
> {
> // Do not apply previous time-step mesh compression flux
> // if the mesh topology changed
> if (mesh.topoChanging())
> {
> talphaPhi1Corr0.clear();
> }
>
> gh = (g & mesh.C()) - ghRef;
> ghf = (g & mesh.Cf()) - ghRef;
>
> MRF.update();
>
> if (correctPhi)
> {
> // Calculate absolute flux
> // from the mapped surface velocity
> phi = mesh.Sf() & Uf();
>
> #include "correctPhi.H"
>
> // Make the flux relative to the mesh motion
> fvc::makeRelative(phi, U);
>
> mixture.correct();
> }
>
> if (checkMeshCourantNo)
> {
> #include "meshCourantNo.H"
> }
> }
> }
>
156a153,157
>
> if (pimple.frozenFlow())
> {
> continue;
> }
とりあえずチュートリアル実行(multiphase/MPPICInterFoam/twoPhasePachuka)
円筒柱形状の計算領域の6割ぐらい水がたまっていて、底面中心の穴(inlet)から、空気(alpha=0)を吹き出すとともに、サイズ分布を持ったパーセルを投入するというものです。
Pachuka Valveとかで製品があるようなので、人の名前でしょうか。エアレーションなどに用いられるようです。
$ cd $FOAM_TURORIALS/multiphase/MPPICInterFoam/twoPhasePachuka
$ cat ./Allrun
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions
# create the underlying block mesh
m4 system/pachuka.m4 > system/blockMeshDict
runApplication blockMesh
\cp 0/alpha.water.orig 0/alpha.water
# create faceSet for burner inlet and faceZone for coupled wall
runApplication topoSet
# create burner inlet
runApplication createPatch -overwrite
# Set alpha.water
runApplication setFields
# Decompose mesh
decomposePar > log.decomposePar 2>&1
# Run
runParallel $(getApplication)
# Reconstruct case
runApplication reconstructPar
$ ./Allrun
constant/kinematicCloudProperties
パーセルの設定となるのが、このファイルですが、設定可能項目が多く、しかも粒子系の専門家でないとわからないモデル名もあって難解です。ネットにも断片的なものしかないので、少し調べてみることにしました。
今回はv1812(or -plus)のものなので、Foundation版やバージョン違いだと、設定可能項目が異なる可能性があります。
下のほうにあるdampingModelとか、isotropyModelとかの指定はどうもよくわかりませんでした。粒子系の物理をやっている人であればわかるのかもしれません。
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Version: v1812
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object kinematicCloudProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// パーセル時間進行計算の設定
solution
{
active true; //Activate or deactivate the particle cloud (yes)
coupled true; //Enable or disable phase coupling (true)
transient yes; //Transient or steady-state solution (yes)
cellValueSourceCorrection no; //Enable/disable correction of momentum transferred to the Eulerian phase (no)
maxCo 0.3; // Co number
interpolationSchemes // Interpolation/integration schemes for the Lagrangian Particle Tracking
{
rho cells;
U cellPoint;
mu cell;
gradAlpha cellPoint;
}
// lagrangian averaging procedure. basic/dual/moment
// basic: cell-volume based average (constant in cell)
// dual: Point values are summed using the tetrahedral decomposition of the computational cells.
// Summation is done in the cells, and also in the terahedrons surrounding each point
// moment: Point values and moments from the cell centroid are summed over
// computational cells. A linear function is generated which has the same
// integrated moment as that of the point data.
averagingMethod dual;
integrationSchemes // time integration scheme
{
U Euler;
}
sourceTerms // source term
{
schemes
{
U semiImplicit 1.0; //explicit or semiImplicit, relaxCoeff for each of the fields
}
}
}
// パーティクルの物理量
constantProperties // Properties of the particles (discrete phase)
{
rho0 1050; // density of particle
alphaMax 0.99; // maximum dispersed phase-fraction (e.g. packing limit)
}
subModels // Specify applicable forces on the particles
{
// パーティクルにかかる力(浮力、抗力、重力など)；必要なものを複数記述する
// MPPICInterFoamではDrag、gravity, interface(気液界面で受ける力）
particleForces
{
/* OF-plus, can be specified multi-forces, for OF-6 different models are included
Drag models: nonSphereDrag
WenYuDrag
ErgunWenYuDrag
PlessisMasliyahDrag
sphereDrag
Lift models: SaffmanMeiLiftForce
TomiyamaLift
SRF ... particle SRF reference frame force (SRF計算する際に必要か？)
gravity ... 重力による沈降など
interface .... Vector force apply to particles to avoid the crossing of particles from
one phase to the other, The magnitude is calculated as C*mass*grad(alpha)
nonInertialFrame ... particle non-inertial reference frame force
paramagnetic ... particle paramagnetic (magnetic field) force
pressureGradient ... particle pressure gradient force
virtualMass .... particle virtual mass force (https://en.wikipedia.org/wiki/Added_mass)
*/
WenYuDrag
{
alphac alphac; // Continuous phase fraction
}
gravity;
interface
{
C -10;
alpha alpha.water;
}
}
// パーティクルのInjectionモデル選択
injectionModels
{
model1 // インジェクション位置の名前
{
/*
type
選択肢：patchFlowRateInjection; patchInjection;
cellZoneInjection; coneInjection, coneNozzleInjection,
fieldActivatedInjection, inflationInjection, injectionModel,
kinematicLookupTableInjection, manualInjection, noInjection
manualInjection can be used for instant injections.
PatchInjection allows time-dependent injection.
*/
// Injectionのタイプにより指定すべきパラメータが異なる.例えばpatchInjectionであれば、↓
type patchInjection; //patchFlowRateInjection;
patch inlet;
SOI 0; // start time of injection, s
duration 10; // duration time, s
parcelsPerSecond 100; // injected parcels per second
U0 (0 0.5 0); // initial parcel/particle velocity
parcelBasisType fixed; // fixed; number;
massTotal 1; // total mass, kg
flowRateProfile constant 1; //
nParticle 1; // number of particle per percel
meanParticleDiameter 0.005;
sizeDistribution // Parcell size variation(9 types:RosinRammler,binned,exponential,fixedValue,general,massRosinRammler,multiNormal,normal,uniform
)
{
type fixedValue; // 固定分布
fixedValueDistribution
{
value 0.005; //m
}
/*
type normal; // 正規分布
normalDistribution
{
expectation 100e-6;
variance 25e-6;
minValue 20e-6;
maxValue 180e-6;
}
type RosinRammler; // RosinRammler分布
RosinRammlerDistribution
{
minValue 200e-6;
maxValue 300e-6;
d 250e-6;
n 3;
}
*/
}
}
}
// パーセルの分散モデル（ソルバーによって実装の有無がある。MPPIC系はnoneのみ）
//Turbulent dispersion models (Discrete Random Walk, Gradient Dispersion)
dispersionModel none;
// パーセルと境界パッチのとの相互作用の指定
// パッチ毎にしていするならlocalInteractionを選択する
patchInteractionModel localInteraction; //standardWallInteraction;
localInteractionCoeffs
{
patches
(
walls
{
type rebound; //rebound（反射）, stick（粘着） or escape（脱出）
e 0.97; // eは、入射速度垂直成分に対する反発係数（1.0完全弾性）
mu 0.09; // muは、入射速度接線成分に対する減衰係数（0.0完全弾性）
}
inlet
{
type rebound;
e 0.97;
mu 0.09;
}
"outlet.*"
{
type escape;
}
);
}
StandardWallInteractionCoeffs
{
type rebound;
e 0.97;
mu 0.09;
}
// パーセルの熱移動モデル（MPPICInterFoamは熱移動は実装されていないのでdummy）
heatTransferModel none;
// Surface film model for dripping and film interaction （MPPICInterFoamは熱移動は実装されていないのでdummy）
surfaceFilmModel none; //absorb, bounce, splash
// パーセルの充填（packing）モデル（セルに対する充填率。パーセルが球と仮定すると最密充填率0.68）
packingModel implicit; //explicit;
explicitCoeffs
{
particleStressModel
{
type HarrisCrighton;
alphaPacked 0.6;
pSolid 10.0;
beta 2.0;
eps 1.0e-7;
}
correctionLimitingMethod
{
type absolute;
e 0.9;
}
}
implicitCoeffs
{
alphaMin 0.001;
rhoMin 1.0;
applyLimiting false;
applyGravity false;
particleStressModel
{
type HarrisCrighton;
alphaPacked 0.6;
pSolid 5.0;
beta 2.0;
eps 1.0e-2;
}
}
// ?
dampingModel relaxation;
relaxationCoeffs
{
timeScaleModel
{
type nonEquilibrium;
alphaPacked 0.7;
e 0.8;
}
}
isotropyModel stochastic;
stochasticCoeffs
{
timeScaleModel
{
type isotropic;
alphaPacked 0.6;
e 0.9;
}
}
stochasticCollisionModel none;
radiation off;
}
cloudFunctions
{
}
// ************************************************************************* //