BUG: compressibleInterFoam family: Corrected transonic option

Resolves bug-report https://bugs.openfoam.org/view.php?id=2785

ENH: compressibleInterFoam family: merged two-phase momentum stress modelling from compressibleInterPhaseTransportFoam

The new momentum stress model selector class
compressibleInterPhaseTransportModel is now used to select between the options:

Description
    Transport model selection class for the compressibleInterFoam family of
    solvers.

    By default the standard mixture transport modelling approach is used in
    which a single momentum stress model (laminar, non-Newtonian, LES or RAS) is
    constructed for the mixture.  However if the \c simulationType in
    constant/turbulenceProperties is set to \c twoPhaseTransport the alternative
    Euler-Euler two-phase transport modelling approach is used in which separate
    stress models (laminar, non-Newtonian, LES or RAS) are instantiated for each
    of the two phases allowing for different modeling for the phases.

Mixture and two-phase momentum stress modelling is now supported in
compressibleInterFoam, compressibleInterDyMFoam and compressibleInterFilmFoam.
The prototype compressibleInterPhaseTransportFoam solver is no longer needed and
has been removed.

applications/solvers/multiphase/VoF/createAlphaFluxes.H

@@ -10,6 +10,11 @@ IOobject alphaPhi10Header
 const bool alphaRestart =
     alphaPhi10Header.typeHeaderOk<surfaceScalarField>(true);
 
+if (alphaRestart)
+{
+    Info << "Restarting alpha" << endl;
+}
+
 // MULES flux from previous time-step
 surfaceScalarField alphaPhi10
 (

applications/solvers/multiphase/compressibleInterFoam/Allwclean

@@ -3,6 +3,8 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 wclean libso twoPhaseMixtureThermo
 wclean libso surfaceTensionModels
+wclean libso VoFphaseCompressibleTurbulenceModels
+
 wclean
 wclean compressibleInterDyMFoam
 wclean compressibleInterFilmFoam

applications/solvers/multiphase/compressibleInterFoam/Allwmake

@@ -6,10 +6,10 @@ cd ${0%/*} || exit 1                        # Run from this directory
 
 wmake $targetType twoPhaseMixtureThermo
 wmake $targetType surfaceTensionModels
+wmake $targetType VoFphaseCompressibleTurbulenceModels
 
 wmake $targetType
 wmake $targetType compressibleInterDyMFoam
 wmake $targetType compressibleInterFilmFoam
-compressibleInterPhaseTransportFoam/Allwmake $targetType $*
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/compressibleInterFoam/Make/options

@@ -8,6 +8,8 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -IVoFphaseCompressibleTurbulenceModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -22,6 +24,7 @@ EXE_LIBS = \
     -linterfaceProperties \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
+    -lVoFphaseCompressibleTurbulenceModels \
     -lfiniteVolume \
     -lfvOptions \
     -lmeshTools

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -3,7 +3,7 @@
     (
         fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
       - fvm::Sp(contErr, T)
-      - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+      - fvm::laplacian(turbulence.alphaEff(), T)
       + (
             divU*p
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -3,7 +3,7 @@
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
       - fvm::Sp(contErr, U)
       + MRF.DDt(rho, U)
-      + turbulence->divDevRhoReff(U)
+      + turbulence.divDevRhoReff(U)
      ==
         fvOptions(rho, U)
     );

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels/Make/files

@@ -1,3 +1,4 @@
 VoFphaseCompressibleTurbulenceModels.C
+compressibleInterPhaseTransportModel.C
 
 LIB = $(FOAM_LIBBIN)/libVoFphaseCompressibleTurbulenceModels

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I../twoPhaseMixtureThermo \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
@@ -9,12 +11,15 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
+    -ltwoPhaseMixtureThermo \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \
+    -ltwoPhaseMixture \
+    -ltwoPhaseProperties \
+    -linterfaceProperties \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lfvOptions \
     -lmeshTools

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels/compressibleInterPhaseTransportModel.C

@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleInterPhaseTransportModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const surfaceScalarField& rhoPhi,
+    const surfaceScalarField& alphaPhi10,
+    const twoPhaseMixtureThermo& mixture
+)
+:
+    twoPhaseTransport_(false),
+    mixture_(mixture),
+    phi_(phi),
+    alphaPhi10_(alphaPhi10)
+{
+    {
+        IOdictionary turbulenceProperties
+        (
+            IOobject
+            (
+                turbulenceModel::propertiesName,
+                U.time().constant(),
+                U.db(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        word simulationType
+        (
+            turbulenceProperties.lookup("simulationType")
+        );
+
+        if (simulationType == "twoPhaseTransport")
+        {
+            twoPhaseTransport_ = true;
+        }
+    }
+
+    if (twoPhaseTransport_)
+    {
+        const volScalarField& alpha1(mixture_.alpha1());
+        const volScalarField& alpha2(mixture_.alpha2());
+
+        const volScalarField& rho1 = mixture_.thermo1().rho();
+        const volScalarField& rho2 = mixture_.thermo2().rho();
+
+        alphaRhoPhi1_ =
+        (
+            new surfaceScalarField
+            (
+                IOobject::groupName("alphaRhoPhi", alpha1.group()),
+                fvc::interpolate(rho1)*alphaPhi10_
+            )
+        );
+
+        alphaRhoPhi2_ =
+        (
+            new surfaceScalarField
+            (
+                IOobject::groupName("alphaRhoPhi", alpha2.group()),
+                fvc::interpolate(rho2)*(phi_ - alphaPhi10_)
+            )
+        );
+
+        turbulence1_ =
+        (
+            PhaseCompressibleTurbulenceModel<fluidThermo>::New
+            (
+                alpha1,
+                rho1,
+                U,
+                alphaRhoPhi1_(),
+                phi,
+                mixture.thermo1()
+            )
+        );
+
+        turbulence2_ =
+        (
+            PhaseCompressibleTurbulenceModel<fluidThermo>::New
+            (
+                alpha2,
+                rho2,
+                U,
+                alphaRhoPhi2_(),
+                phi,
+                mixture.thermo2()
+            )
+        );
+    }
+    else
+    {
+        turbulence_ = compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            rhoPhi,
+            mixture
+        );
+
+        turbulence_->validate();
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::compressibleInterPhaseTransportModel::alphaEff() const
+{
+    if (twoPhaseTransport_)
+    {
+        return
+            mixture_.alpha1()*mixture_.thermo1().alphaEff(turbulence1_->mut())
+          + mixture_.alpha2()*mixture_.thermo2().alphaEff(turbulence2_->mut());
+    }
+    else
+    {
+        return turbulence_->mut();
+    }
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::compressibleInterPhaseTransportModel::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    if (twoPhaseTransport_)
+    {
+        return
+            turbulence1_->divDevRhoReff(U)
+          + turbulence2_->divDevRhoReff(U);
+    }
+    else
+    {
+        return turbulence_->divDevRhoReff(U);
+    }
+}
+
+
+void Foam::compressibleInterPhaseTransportModel::correctPhasePhi()
+{
+    if (twoPhaseTransport_)
+    {
+        const volScalarField& rho1 = mixture_.thermo1().rho();
+        const volScalarField& rho2 = mixture_.thermo2().rho();
+
+        alphaRhoPhi1_.ref() = fvc::interpolate(rho1)*alphaPhi10_;
+        alphaRhoPhi2_.ref() = fvc::interpolate(rho2)*(phi_ - alphaPhi10_);
+    }
+}
+
+
+void Foam::compressibleInterPhaseTransportModel::correct()
+{
+    if (twoPhaseTransport_)
+    {
+        turbulence1_->correct();
+        turbulence2_->correct();
+    }
+    else
+    {
+        turbulence_->correct();
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels/compressibleInterPhaseTransportModel.H

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressibleInterPhaseTransportModel
+
+Description
+    Transport model selection class for the compressibleInterFoam family of
+    solvers.
+
+    By default the standard mixture transport modelling approach is used in
+    which a single momentum stress model (laminar, non-Newtonian, LES or RAS) is
+    constructed for the mixture.  However if the \c simulationType in
+    constant/turbulenceProperties is set to \c twoPhaseTransport the alternative
+    Euler-Euler two-phase transport modelling approach is used in which separate
+    stress models (laminar, non-Newtonian, LES or RAS) are instantiated for each
+    of the two phases allowing for different modeling for the phases.
+
+SourceFiles
+    compressibleInterPhaseTransportModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleInterPhaseTransportModel_H
+#define compressibleInterPhaseTransportModel_H
+
+#include "twoPhaseMixture.H"
+#include "twoPhaseMixtureThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "VoFphaseCompressibleTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class compressibleInterPhaseTransportModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class compressibleInterPhaseTransportModel
+{
+    // Private data
+
+        //- Switch to select two-phase or mixture transport modelling
+        Switch twoPhaseTransport_;
+
+        //- Two-phase mixture
+        const twoPhaseMixtureThermo& mixture_;
+
+        //- Mixture volumetric flux
+        const surfaceScalarField& phi_;
+
+        //- Phase volumetric flux
+        const surfaceScalarField& alphaPhi10_;
+
+        //- Phase-1 mass-flux (constructed for two-phase transport)
+        tmp<surfaceScalarField> alphaRhoPhi1_;
+
+        //- Phase-2 mass-flux (constructed for two-phase transport)
+        tmp<surfaceScalarField> alphaRhoPhi2_;
+
+        //- Mixture transport model (constructed for mixture transport)
+        autoPtr<compressible::turbulenceModel> turbulence_;
+
+        //- Phase-1 transport model (constructed for two-phase transport)
+        autoPtr<PhaseCompressibleTurbulenceModel<fluidThermo>> turbulence1_;
+
+        //- Phase-2 transport model (constructed for two-phase transport)
+        autoPtr<PhaseCompressibleTurbulenceModel<fluidThermo>> turbulence2_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        compressibleInterPhaseTransportModel
+        (
+            const compressibleInterPhaseTransportModel&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const compressibleInterPhaseTransportModel&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        compressibleInterPhaseTransportModel
+        (
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const surfaceScalarField& rhoPhi,
+            const surfaceScalarField& alphaPhi10,
+            const twoPhaseMixtureThermo& mixture
+        );
+
+
+    // Member Functions
+
+        //- Return the effective temperature transport coefficient
+        tmp<volScalarField> alphaEff() const;
+
+        //- Return the effective momentum stress divergence
+        tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+
+        //- Correct the phase mass-fluxes
+        // (required for the two-phase transport option)
+        void correctPhasePhi();
+
+        //- Correct the phase or mixture transport models
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options

@@ -3,12 +3,14 @@ EXE_INC = \
     -I.. \
     -I../../VoF \
     -I../twoPhaseMixtureThermo \
+    -I../VoFphaseCompressibleTurbulenceModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
@@ -25,6 +27,7 @@ EXE_LIBS = \
     -linterfaceProperties \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
+    -lVoFphaseCompressibleTurbulenceModels \
     -ldynamicMesh \
     -lmeshTools \
     -ldynamicFvMesh \

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -44,9 +44,7 @@ Description
 #include "localEulerDdtScheme.H"
 #include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
-#include "twoPhaseMixture.H"
-#include "twoPhaseMixtureThermo.H"
-#include "turbulentFluidThermoModel.H"
+#include "compressibleInterPhaseTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
@@ -64,7 +62,6 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createControl.H"
     #include "createFields.H"
-    #include "createAlphaFluxes.H"
     #include "createUf.H"
     #include "createControls.H"
     #include "CourantNo.H"
@@ -75,8 +72,6 @@ int main(int argc, char *argv[])
     const volScalarField& psi1 = mixture.thermo1().psi();
     const volScalarField& psi2 = mixture.thermo2().psi();
 
-    turbulence->validate();
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
 
@@ -149,6 +144,8 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
+            turbulence.correctPhasePhi();
+
             #include "UEqn.H"
             #include "TEqn.H"
 
@@ -160,7 +157,7 @@ int main(int argc, char *argv[])
 
             if (pimple.turbCorr())
             {
-                turbulence->correct();
+                turbulence.correct();
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -31,27 +31,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
         surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
 
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
-          + correction
-            (
-                psi1*fvm::ddt(p_rgh)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
-          + correction
-            (
-                psi2*fvm::ddt(p_rgh)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
-            );
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options

@@ -3,12 +3,14 @@ EXE_INC = \
     -I.. \
     -I../../VoF \
     -I../twoPhaseMixtureThermo \
+    -I../VoFphaseCompressibleTurbulenceModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
@@ -34,6 +36,7 @@ EXE_LIBS = \
     -linterfaceProperties \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
+    -lVoFphaseCompressibleTurbulenceModels \
     -lSLGThermo \
     -lsurfaceFilmModels \
     -lsurfaceFilmDerivedFvPatchFields \

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/TEqn.H

@@ -3,7 +3,7 @@
     (
         fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
       - fvm::Sp(contErr, T)
-      - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+      - fvm::laplacian(turbulence.alphaEff(), T)
       + (
             (
                 fvc::div(fvc::absolute(phi, U), p)

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C

@@ -41,10 +41,8 @@ Description
 #include "localEulerDdtScheme.H"
 #include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
-#include "rhoThermo.H"
-#include "twoPhaseMixture.H"
-#include "twoPhaseMixtureThermo.H"
-#include "turbulentFluidThermoModel.H"
+#include "compressibleInterPhaseTransportModel.H"
+#include "pimpleControl.H"
 #include "SLGThermo.H"
 #include "surfaceFilmModel.H"
 #include "pimpleControl.H"
@@ -63,7 +61,6 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createFields.H"
-    #include "createAlphaFluxes.H"
     #include "createSurfaceFilmModel.H"
 
     volScalarField& p = mixture.p();
@@ -73,8 +70,6 @@ int main(int argc, char *argv[])
 
     regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
 
-    turbulence->validate();
-
     if (!LTS)
     {
         #include "readTimeControls.H"
@@ -113,6 +108,8 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
+            turbulence.correctPhasePhi();
+
             volScalarField::Internal Srho(surfaceFilm.Srho());
             contErr -= posPart(Srho);
 
@@ -127,7 +124,7 @@ int main(int argc, char *argv[])
 
             if (pimple.turbCorr())
             {
-                turbulence->correct();
+                turbulence.correct();
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H

@@ -28,27 +28,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
         surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
 
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
-          + correction
-            (
-                psi1*fvm::ddt(p_rgh)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
-          + correction
-            (
-                psi2*fvm::ddt(p_rgh)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
-            );
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -34,7 +34,10 @@ Description
     The momentum and other fluid properties are of the "mixture" and a single
     momentum equation is solved.
 
-    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+    Either mixture or two-phase transport modelling may be selected.  In the
+    mixture approach a single laminar, RAS or LES model is selected to model the
+    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
+    RAS or LES selected models are selected for each of the phases.
 
 \*---------------------------------------------------------------------------*/
 
@@ -44,10 +47,7 @@ Description
 #include "localEulerDdtScheme.H"
 #include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
-#include "rhoThermo.H"
-#include "twoPhaseMixture.H"
-#include "twoPhaseMixtureThermo.H"
-#include "turbulentFluidThermoModel.H"
+#include "compressibleInterPhaseTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "fvcSmooth.H"
@@ -64,15 +64,12 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createFields.H"
-    #include "createAlphaFluxes.H"
 
     volScalarField& p = mixture.p();
     volScalarField& T = mixture.T();
     const volScalarField& psi1 = mixture.thermo1().psi();
     const volScalarField& psi2 = mixture.thermo2().psi();
 
-    turbulence->validate();
-
     if (!LTS)
     {
         #include "readTimeControls.H"
@@ -109,6 +106,8 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
+            turbulence.correctPhasePhi();
+
             #include "UEqn.H"
             volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
@@ -121,7 +120,7 @@ int main(int argc, char *argv[])
 
             if (pimple.turbCorr())
             {
-                turbulence->correct();
+                turbulence.correct();
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/Allwclean

@@ -1,7 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-wclean libso VoFphaseCompressibleTurbulenceModels
-wclean
-
-#------------------------------------------------------------------------------

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/Allwmake

@@ -1,10 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Parse arguments for library compilation
-. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-
-wmake $targetType VoFphaseCompressibleTurbulenceModels
-wmake $targetType
-
-#------------------------------------------------------------------------------

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/Make/files

@@ -1,3 +0,0 @@
-compressibleInterPhaseTransportFoam.C
-
-EXE = $(FOAM_APPBIN)/compressibleInterPhaseTransportFoam

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/Make/options

@@ -1,31 +0,0 @@
-EXE_INC = \
-    -I. \
-    -I$(FOAM_SOLVERS)/multiphase/VoF \
-    -I$(FOAM_SOLVERS)/multiphase/compressibleInterFoam/twoPhaseMixtureThermo \
-    -I$(FOAM_SOLVERS)/multiphase/compressibleInterFoam \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
-    -IVoFphaseCompressibleTurbulenceModels/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -ltwoPhaseMixtureThermo \
-    -ltwoPhaseSurfaceTension \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -ltwoPhaseMixture \
-    -ltwoPhaseProperties \
-    -linterfaceProperties \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lVoFphaseCompressibleTurbulenceModels \
-    -lfiniteVolume \
-    -lfvOptions \
-    -lmeshTools

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/TEqn.H

@@ -1,37 +0,0 @@
-{
-    fvScalarMatrix TEqn
-    (
-        fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
-      - fvm::Sp(contErr, T)
-      - fvm::laplacian
-        (
-            mixture.alphaEff
-            (
-               alpha1*turbulence1->mut()
-             + alpha2*turbulence2->mut()
-            ),
-            T
-        )
-      + (
-            fvc::div(fvc::absolute(phi, U), p)
-          + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
-        )
-       *(
-           alpha1/mixture.thermo1().Cv()
-         + alpha2/mixture.thermo2().Cv()
-        )
-     ==
-        fvOptions(rho, T)
-    );
-
-    TEqn.relax();
-
-    fvOptions.constrain(TEqn);
-
-    TEqn.solve();
-
-    fvOptions.correct(T);
-
-    mixture.correctThermo();
-    mixture.correct();
-}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/UEqn.H

@@ -1,35 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      - fvm::Sp(contErr, U)
-      + MRF.DDt(rho, U)
-      + turbulence1->divDevRhoReff(U)
-      + turbulence2->divDevRhoReff(U)
-     ==
-        fvOptions(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
-
-        fvOptions.correct(U);
-
-        K = 0.5*magSqr(U);
-    }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/compressibleInterPhaseTransportFoam.C

@@ -1,180 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    compressibleInterPhaseTransportFoam
-
-Description
-    Solver for 2 compressible, non-isothermal immiscible fluids using a VOF
-    (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.
-
-    The fluid stress modelling is generic Euler-Euler two-phase in which
-    separate laminar, RAS or LES are selected for each of the phases.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "CMULES.H"
-#include "EulerDdtScheme.H"
-#include "localEulerDdtScheme.H"
-#include "CrankNicolsonDdtScheme.H"
-#include "subCycle.H"
-#include "rhoThermo.H"
-#include "twoPhaseMixture.H"
-#include "twoPhaseMixtureThermo.H"
-#include "turbulentFluidThermoModel.H"
-#include "pimpleControl.H"
-#include "fvOptions.H"
-#include "fvcSmooth.H"
-#include "VoFphaseCompressibleTurbulenceModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
-    #include "createFields.H"
-    #include "createAlphaFluxes.H"
-
-    volScalarField& p = mixture.p();
-    volScalarField& T = mixture.T();
-    const volScalarField& psi1 = mixture.thermo1().psi();
-    const volScalarField& psi2 = mixture.thermo2().psi();
-
-
-    surfaceScalarField alphaRhoPhi1
-    (
-        IOobject::groupName("alphaRhoPhi", alpha1.group()),
-        fvc::interpolate(rho1)*alphaPhi10
-    );
-
-    autoPtr<PhaseCompressibleTurbulenceModel<fluidThermo>> turbulence1
-    (
-        PhaseCompressibleTurbulenceModel<fluidThermo>::New
-        (
-            alpha1,
-            rho1,
-            U,
-            alphaRhoPhi1,
-            phi,
-            mixture.thermo1()
-        )
-    );
-
-    surfaceScalarField alphaRhoPhi2
-    (
-        IOobject::groupName("alphaRhoPhi", alpha2.group()),
-        fvc::interpolate(rho2)*(phi - alphaPhi10)
-    );
-
-    autoPtr<PhaseCompressibleTurbulenceModel<fluidThermo>> turbulence2
-    (
-        PhaseCompressibleTurbulenceModel<fluidThermo>::New
-        (
-            alpha2,
-            rho2,
-            U,
-            alphaRhoPhi2,
-            phi,
-            mixture.thermo2()
-        )
-    );
-
-    if (!LTS)
-    {
-        #include "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readTimeControls.H"
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "CourantNo.H"
-            #include "alphaCourantNo.H"
-            #include "setDeltaT.H"
-        }
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "alphaControls.H"
-            #include "compressibleAlphaEqnSubCycle.H"
-
-            alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi10;
-            alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi10);
-
-            #include "UEqn.H"
-            #include "TEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence1->correct();
-                turbulence2->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = "
-            << runTime.elapsedCpuTime()
-            << " s\n\n" << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterPhaseTransportFoam/createFields.H

@@ -1,96 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading field p_rgh\n" << endl;
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
-twoPhaseMixtureThermo mixture(U, phi);
-
-volScalarField& alpha1(mixture.alpha1());
-volScalarField& alpha2(mixture.alpha2());
-
-Info<< "Reading thermophysical properties\n" << endl;
-
-const volScalarField& rho1 = mixture.thermo1().rho();
-const volScalarField& rho2 = mixture.thermo2().rho();
-
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    alpha1*rho1 + alpha2*rho2
-);
-
-
-dimensionedScalar pMin
-(
-    "pMin",
-    dimPressure,
-    mixture
-);
-
-mesh.setFluxRequired(p_rgh.name());
-mesh.setFluxRequired(alpha1.name());
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-// Mass flux
-// Initialisation does not matter because rhoPhi is reset after the
-// alpha1 solution before it is used in the U equation.
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    fvc::interpolate(rho)*phi
-);
-
-volScalarField dgdt(alpha1*fvc::div(phi));
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-#include "createMRF.H"
-#include "createFvOptions.H"

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -89,10 +89,17 @@ surfaceScalarField rhoPhi
 
 volScalarField dgdt(alpha1*fvc::div(phi));
 
+#include "createAlphaFluxes.H"
+
 // Construct compressible turbulence model
-autoPtr<compressible::turbulenceModel> turbulence
+compressibleInterPhaseTransportModel turbulence
 (
-    compressible::turbulenceModel::New(rho, U, rhoPhi, mixture)
+    rho,
+    U,
+    phi,
+    rhoPhi,
+    alphaPhi10,
+    mixture
 );
 
 #include "createK.H"

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -28,27 +28,43 @@
 
     if (pimple.transonic())
     {
+        #include "rhofs.H"
+
         surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
         surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
 
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
-          + correction
-            (
-                psi1*fvm::ddt(p_rgh)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1/rho1)
+               *correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                )
             );
-        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
         p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
-          + correction
-            (
-                psi2*fvm::ddt(p_rgh)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
-            );
-        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2/rho2)
+              *correction
+               (
+                   psi2*fvm::ddt(p_rgh)
+                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+               )
+           );
         p_rghEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -326,10 +326,6 @@ while (pimple.correct())
                     ).ptr()
                 );
 
-                deleteDemandDrivenData
-                (
-                    pEqnComps[phasei].faceFluxCorrectionPtr()
-                );
                 pEqnComps[phasei].relax();
             }
             else

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -255,7 +255,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
         else
@@ -297,7 +296,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
             pEqnComp2.ref().relax();
         }
         else

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -241,7 +241,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
         else
@@ -283,7 +282,6 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
             pEqnComp2.ref().relax();
         }
         else

applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H

@@ -262,7 +262,6 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
         pEqnComp1.ref().relax();
 
         pEqnComp2 =
@@ -278,7 +277,6 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
         pEqnComp2.ref().relax();
     }
     else

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H

@@ -243,7 +243,6 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
         pEqnComp1.ref().relax();
 
         pEqnComp2 =
@@ -259,7 +258,6 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
         pEqnComp2.ref().relax();
     }
     else