From 5d3d40392f73acede477a59ea4a66c585665e732 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 25 Jun 2015 22:29:08 +0100
Subject: [PATCH] reactingTwoPhaseEulerFoam: Add fvOption handling to the
 continuity error correction in MovingPhaseModel<BasePhaseModel>::correct()

---
 .../multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H | 7 -------
 .../solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H | 2 +-
 .../multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H        | 2 +-
 .../phaseModel/MovingPhaseModel/MovingPhaseModel.C         | 5 ++++-
 .../phaseSystems/phaseSystem/phaseSystem.H                 | 4 ++--
 .../phaseSystems/phaseSystem/phaseSystemI.H                | 2 +-
 .../reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C  | 2 --
 7 files changed, 9 insertions(+), 15 deletions(-)
 delete mode 100644 applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H

diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H
deleted file mode 100644
index 6aeb5d024bf..00000000000
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/correctContErrs.H
+++ /dev/null
@@ -1,7 +0,0 @@
-phase1.continuityError() =
-    fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
-  - (fvOptions(alpha1, rho1)&rho1);
-
-phase2.continuityError() =
-    fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
-  - (fvOptions(alpha2, rho2)&rho2);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 338b0d88b67..b9029e4190f 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -89,7 +89,7 @@ tmp<surfaceScalarField> phiF2;
 while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
-    #include "correctContErrs.H"
+    fluid.correct();
 
     // Correct fixed-flux BCs to be consistent with the velocity BCs
     MRF.correctBoundaryFlux(U1, phi1);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 8736d0d5a3f..8bd69f606c3 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -93,7 +93,7 @@ tmp<surfaceScalarField> Ff2;
 while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
-    #include "correctContErrs.H"
+    fluid.correct();
 
     surfaceScalarField rhof1(fvc::interpolate(rho1));
     surfaceScalarField rhof2(fvc::interpolate(rho2));
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 8adb8674016..3d2e2750806 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -217,8 +217,11 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correct()
 
     this->fluid().MRF().correctBoundaryVelocity(U_);
 
+    volScalarField& rho = this->thermo().rho();
+
     continuityError_ =
-        fvc::ddt(*this, this->thermo().rho()) + fvc::div(alphaRhoPhi_);
+        fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_)
+      - (this->fluid().fvOptions()(*this, rho)&rho);
 }
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index e12f52df926..691a663ffbe 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -179,7 +179,7 @@ protected:
         IOMRFZoneList MRF_;
 
         //- Optional FV-options
-        fv::IOoptionList fvOptions_;
+        mutable fv::IOoptionList fvOptions_;
 
         //- Blending methods
         blendingMethodTable blendingMethods_;
@@ -371,7 +371,7 @@ public:
             inline const IOMRFZoneList& MRF() const;
 
             //- Optional FV-options
-            inline fv::IOoptionList& fvOptions();
+            inline fv::IOoptionList& fvOptions() const;
 
             //- Access a sub model between a phase pair
             template <class modelType>
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 0613ba47acb..3d7e837533a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
@@ -61,7 +61,7 @@ inline const Foam::IOMRFZoneList& Foam::phaseSystem::MRF() const
 }
 
 
-inline Foam::fv::IOoptionList& Foam::phaseSystem::fvOptions()
+inline Foam::fv::IOoptionList& Foam::phaseSystem::fvOptions() const
 {
     return fvOptions_;
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index 087434ba5b1..f74e3e93f08 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -94,8 +94,6 @@ int main(int argc, char *argv[])
             fluid.solve();
             fluid.correct();
 
-            #include "correctContErrs.H"
-
             #include "YEqns.H"
 
             if (faceMomentum)
-- 
GitLab