twoPhaseEulerFoam: Add MRF support and remove MRFTwoPhaseEulerFoam

applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean

@@ -6,6 +6,5 @@ wclean libso phaseModel
 wclean libso interfacialModels
 wclean libso kineticTheoryModels
 wclean
-wclean MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake

@@ -6,6 +6,5 @@ wmake libso phaseModel
 wmake libso interfacialModels
 wmake libso kineticTheoryModels
 wmake
-wmake MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C

@@ -1,119 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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
-    twoPhaseEulerFoam
-
-Description
-    Solver for a system of 2 incompressible fluid phases with one phase
-    dispersed, e.g. gas bubbles in a liquid or solid particles in a gas.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "nearWallDist.H"
-#include "wallFvPatch.H"
-#include "Switch.H"
-
-#include "IFstream.H"
-#include "OFstream.H"
-
-#include "dragModel.H"
-#include "phaseModel.H"
-#include "kineticTheoryModel.H"
-
-#include "pimpleControl.H"
-#include "MRFZones.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
-    #include "createFields.H"
-    #include "readPPProperties.H"
-    #include "initContinuityErrs.H"
-    #include "createMRFZones.H"
-    #include "readTimeControls.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    pimpleControl pimple(mesh);
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readTwoPhaseEulerFoamControls.H"
-        #include "CourantNos.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "alphaEqn.H"
-            #include "liftDragCoeffs.H"
-            #include "UEqns.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-
-                if (correctAlpha && !pimple.finalIter())
-                {
-                    #include "alphaEqn.H"
-                }
-            }
-
-            #include "DDtU.H"
-
-            if (pimple.turbCorr())
-            {
-                #include "kEpsilon.H"
-            }
-        }
-
-        #include "write.H"
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files

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

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options

@@ -1,17 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -I../../bubbleFoam \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -IturbulenceModel \
-    -I../kineticTheoryModels/lnInclude \
-    -I../interfacialModels/lnInclude \
-    -I../phaseModel/lnInclude \
-
-EXE_LIBS = \
-    -lEulerianInterfacialModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lincompressibleTransportModels \
-    -lphaseModel \
-    -lkineticTheoryModel

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H

@@ -1,99 +0,0 @@
-fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
-fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
-
-{
-    {
-        volTensorField gradU1T(T(fvc::grad(U1)));
-
-        if (kineticTheory.on())
-        {
-            kineticTheory.solve(gradU1T);
-            nuEff1 = kineticTheory.mu1()/rho1;
-        }
-        else // If not using kinetic theory is using Ct model
-        {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
-        }
-
-        volTensorField Rc1
-        (
-            "Rc1",
-            (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
-        );
-
-        if (kineticTheory.on())
-        {
-            Rc1 -= ((kineticTheory.lambda()/rho1)*tr(gradU1T))*tensor(I);
-        }
-
-        surfaceScalarField phiR1
-        (
-            -fvc::interpolate(nuEff1)*mesh.magSf()*fvc::snGrad(alpha1)
-            /fvc::interpolate(alpha1 + scalar(0.001))
-        );
-
-        U1Eqn =
-        (
-            (scalar(1) + Cvm*rho2*alpha2/rho1)*
-            (
-                fvm::ddt(U1)
-              + fvm::div(phi1, U1, "div(phi1,U1)")
-              - fvm::Sp(fvc::div(phi1), U1)
-            )
-
-          - fvm::laplacian(nuEff1, U1)
-          + fvc::div(Rc1)
-
-          + fvm::div(phiR1, U1, "div(phi1,U1)")
-          - fvm::Sp(fvc::div(phiR1), U1)
-          + (fvc::grad(alpha1)/(fvc::average(alpha1) + scalar(0.001)) & Rc1)
-         ==
-        //  g                          // Buoyancy term transfered to p-equation
-          - fvm::Sp(alpha2/rho1*K, U1)
-        //+ alpha2/rho1*K*U2           // Explicit drag transfered to p-equation
-          - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
-        );
-        mrfZones.addCoriolis(U1Eqn);
-        U1Eqn.relax();
-    }
-
-    {
-        volTensorField gradU2T(T(fvc::grad(U2)));
-        volTensorField Rc2
-        (
-            "Rc2",
-            (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
-        );
-
-        surfaceScalarField phiR2
-        (
-            -fvc::interpolate(nuEff2)*mesh.magSf()*fvc::snGrad(alpha2)
-            /fvc::interpolate(alpha2 + scalar(0.001))
-        );
-
-        U2Eqn =
-        (
-            (scalar(1) + Cvm*rho2*alpha1/rho2)*
-            (
-                fvm::ddt(U2)
-              + fvm::div(phi2, U2, "div(phi2,U2)")
-              - fvm::Sp(fvc::div(phi2), U2)
-            )
-
-          - fvm::laplacian(nuEff2, U2)
-          + fvc::div(Rc2)
-
-          + fvm::div(phiR2, U2, "div(phi2,U2)")
-          - fvm::Sp(fvc::div(phiR2), U2)
-
-          + (fvc::grad(alpha2)/(fvc::average(alpha2) + scalar(0.001)) & Rc2)
-         ==
-        //  g                          // Buoyancy term transfered to p-equation
-          - fvm::Sp(alpha1/rho2*K, U2)
-        //+ alpha1/rho2*K*U1           // Explicit drag transfered to p-equation
-          + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
-        );
-        mrfZones.addCoriolis(U2Eqn);
-        U2Eqn.relax();
-    }
-}

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H

@@ -1,4 +0,0 @@
-    MRFZones mrfZones(mesh);
-    mrfZones.correctBoundaryVelocity(U1);
-    mrfZones.correctBoundaryVelocity(U2);
-    mrfZones.correctBoundaryVelocity(U);

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

@@ -1,118 +0,0 @@
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    rAU1f = fvc::interpolate(rAU1);
-    surfaceScalarField rAU2f(fvc::interpolate(rAU2));
-
-    volVectorField HbyA1("HbyA1", U1);
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2("HbyA2", U2);
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
-    surfaceScalarField phiDrag1
-    (
-        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
-    );
-
-    if (g0.value() > 0.0)
-    {
-        phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
-    }
-
-
-    surfaceScalarField phiDrag2
-    (
-        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
-    );
-
-    // Fix for gravity on outlet boundary.
-    forAll(p.boundaryField(), patchi)
-    {
-        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
-        {
-            phiDrag1.boundaryField()[patchi] = 0.0;
-            phiDrag2.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
-    );
-    mrfZones.relativeFlux(phiHbyA1);
-
-    surfaceScalarField phiHbyA2
-    (
-        "phiHbyA2",
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
-    );
-    mrfZones.relativeFlux(phiHbyA2);
-
-    mrfZones.relativeFlux(phi1.oldTime());
-    mrfZones.relativeFlux(phi1);
-    mrfZones.relativeFlux(phi2.oldTime());
-    mrfZones.relativeFlux(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            surfaceScalarField SfGradp(pEqn.flux()/Dp);
-
-            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            p.relax();
-            SfGradp = pEqn.flux()/Dp;
-
-            U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1);
-            U1.correctBoundaryConditions();
-
-            U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2);
-            U2.correctBoundaryConditions();
-
-            U = alpha1*U1 + alpha2*U2;
-        }
-    }
-}
-
-#include "continuityErrs.H"

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -53,7 +53,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         //+ alpha2/rho1*K*U2           // Explicit drag transfered to p-equation
           - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
         );
-
+        mrfZones.addCoriolis(U1Eqn);
         U1Eqn.relax();
     }
 
@@ -93,7 +93,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         //+ alpha1/rho2*K*U1           // Explicit drag transfered to p-equation
           + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
         );
-
+        mrfZones.addCoriolis(U2Eqn);
         U2Eqn.relax();
     }
 }

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -93,16 +93,22 @@
 
     dimensionedScalar Cvm
     (
+        "Cvm",
+        dimless,
         transportProperties.lookup("Cvm")
     );
 
     dimensionedScalar Cl
     (
+        "Cl",
+        dimless,
         transportProperties.lookup("Cl")
     );
 
     dimensionedScalar Ct
     (
+        "Ct",
+        dimless,
         transportProperties.lookup("Ct")
     );
 

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -14,6 +14,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    mrfZones.absoluteFlux(phi1.oldTime());
+    mrfZones.absoluteFlux(phi1);
+
+    mrfZones.absoluteFlux(phi2.oldTime());
+    mrfZones.absoluteFlux(phi2);
+
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
@@ -29,6 +35,7 @@
         phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
     }
 
+
     surfaceScalarField phiDrag2
     (
         fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
@@ -51,6 +58,7 @@
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
       + phiDrag1
     );
+    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
@@ -59,6 +67,12 @@
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
       + phiDrag2
     );
+    mrfZones.relativeFlux(phiHbyA2);
+
+    mrfZones.relativeFlux(phi1.oldTime());
+    mrfZones.relativeFlux(phi1);
+    mrfZones.relativeFlux(phi2.oldTime());
+    mrfZones.relativeFlux(phi2);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
@@ -83,14 +97,8 @@
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
-            phi1.boundaryField() ==
-                (fvc::interpolate(U1) & mesh.Sf())().boundaryField();
             phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-
-            phi2.boundaryField() ==
-                (fvc::interpolate(U2) & mesh.Sf())().boundaryField();
             phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-
             phi = alpha1f*phi1 + alpha2f*phi2;
 
             p.relax();

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ Description
 #include "kineticTheoryModel.H"
 
 #include "pimpleControl.H"
+#include "MRFZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +57,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "readPPProperties.H"
     #include "initContinuityErrs.H"
+    #include "createMRFZones.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"