MPPICFoamにMRF項を追加してみる①
MPPICFoamで計算したい問題があってケースファイルを用意したところで、MRFが効かないことに気づき悲しい思いをした。
なければ作ればいい、がオープンソースなので改造をやってみる。
ソースコードを取ってくる
cp /opt/OpenFOAM/OpenFOAM-v1806/applications/solvers/lagrangian/DPMFoam/ ./ -r
慣習でとってきたDPMFoamをmyDPMFoamに変更する。
mv -f DPMFoam myDPMFoam
DPMFoamフォルダには派生ソルバも含まれている。それらのソルバもすべてmy~に書き換える。
DPMTurbulenceModelsフォルダは今回操作しないのでいじらない。
~.Cファイルの中にもmyをつける箇所があるのでそこも書き換える。
#define MPPIC
#include "myDPMDyMFoam.C"
後はコンパイルのためにMakeフォルダ下のfilesで指定している名前にもmyをつける。また、コンパイル後の実行ファイルの保存先はFOAM_USER_APPBINにしておく。デフォルトのままだと、元々あるソルバを上書きする可能性がある。
myDPMFoam.C
EXE = $(FOAM_USER_APPBIN)/myDPMFoam
後はインクルードするファイルを記述したMake/optionsを書き換える。例えば、myDPMFoam下のMake/optionsでは同じ階層のDPMTurbulenceModelsを呼び出しているが、さっきmy~に書き換えたので、合わせて書き換える。
EXE_INC = \
-I./myDPMTurbulenceModels/lnInclude \
~
EXE_LIBS = \
~
-lmyDPMTurbulenceModels \
~
最後にAllwmakeを書き換える。DPMTurbulenceModelsとDyMソルバはコンパイルしないことにした。
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
#------------------------------------------------------------------------------
#wmake $targetType DPMTurbulenceModels
wmake $targetType
wmake $targetType myMPPICFoam
#wmake $targetType myDPMDyMFoam
#wmake $targetType myDPMDyMFoam/myMPPICDyMFoam
#------------------------------------------------------------------------------
書き換えが終わったらAllwmakeを実行してコンパイルできるか確認する。エラーが出たら対応する。
まずはソースコードを眺めてみる
myDPMFoamの中身を見ると、pEqn.HとUcEqn.Hがいることがわかる。pimpleFoamのソースコードを見たことがある人は既視感があると思うが、中身も既視感がある。そのため、pimpleFoamのMRFの部分を移植する。移植する際はソルバによってUがUcになっていたりするので、移植先のコードに合わせる。
とりあえずコードを書き換えてコンパイル
pEqn.HとUEqn.Hを以下のように書き換えた。
myDPMFoam/pEqn.H
{
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(HbyA)
+ alphacf*rAUcf*fvc::ddtCorr(Uc, phic)
+ MRF.zeroFilter(fvc::interpolate(rAUc)*fvc::ddtCorr(Uc, phic)) //add
)
);
MRF.makeRelative(phiHbyA); //add
if (p.needReference())
{
adjustPhi(phiHbyA, Uc, p);
}
phiHbyA += phicForces;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, Uc, phiHbyA, rAUcf, MRF); //mod
~
myDPMFoam/UEqn.H
MRF.correctBoundaryVelocity(Uc); //add
fvVectorMatrix UcEqn
(
fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
+ continuousPhaseTurbulence->divDevRhoReff(Uc)
+ MRF.DDt(Uc)
==
(1.0/rhoc)*cloudSU
);
~
myDPMFoam/myDPMDyMFoam/pEqn.H
{
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(HbyA)
+ alphacf*rAUcf*fvc::ddtCorr(Uc, Ucf)
+ MRF.zeroFilter(fvc::interpolate(rAUc)*fvc::ddtCorr(Uc, Ucf)) //add
)
);
MRF.makeRelative(phiHbyA); //add
if (p.needReference())
{
fvc::makeRelative(phiHbyA, Uc);
adjustPhi(phiHbyA, Uc, p);
fvc::makeAbsolute(phiHbyA, Uc);
}
phiHbyA += phicForces;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, Uc, phiHbyA, rAUcf, MRF); //mod
~
myDPMFoam/myDPMFoam.C
~
while (pimple.loop())
{
//add
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
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);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
~
myDPMFoam/createFields.H
~
#include "createMRF.H"
とりあえずコンパイルできた。
検証
検証は変更した点を確認すればいいので、粒子計算よりも、MRFの影響を見る。mixerVessel2Dとかを流用して確認すればいいのではないかと思う。
検証は出来次第追加します。