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とかを流用して確認すればいいのではないかと思う。

検証は出来次第追加します。