VoF solvers: phase-fraction equation and move mesh motion into PIMPLE loop

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -70,7 +70,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -83,14 +83,14 @@ Foam::multiphaseMixtureThermo::multiphaseMixtureThermo
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
     ),
 
     alphas_

applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C

@@ -77,19 +77,26 @@ int main(int argc, char *argv[])
 
         #include "setrDeltaT.H"
 
-        twoPhaseProperties.correct();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        #define LTSSOLVE
-        #include "alphaEqnSubCycle.H"
-        #undef LTSSOLVE
-
-        interface.correct();
-
-        turbulence->correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #define LTSSOLVE
+                #include "alphaEqnSubCycle.H"
+                #undef LTSSOLVE
+
+                interface.correct();
+            }
+
+            turbulence->correct();
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C

@@ -81,15 +81,22 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-        #include "zonePhaseVolumes.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+                #include "zonePhaseVolumes.H"
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interFoam/Make/options

@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -ggdb3 \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -6,9 +6,7 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
-    tmp<surfaceScalarField> tphiAlpha;
+    if (pimple.firstIter() && MULESCorr)
-
-    if (MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -1,5 +1,3 @@
-#include "alphaControls.H"
-
 if (nAlphaSubCycles > 1)
 {
     dimensionedScalar totalDeltaT = runTime.deltaT();

applications/solvers/multiphase/interFoam/createFields.H

@@ -62,13 +62,14 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        rho1*phi
+        mesh,
+        dimensionedScalar("0", dimMass/dimTime, 0)
     );
 
 

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -77,46 +77,56 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        mesh.update();
-
-        if (mesh.changing())
-        {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
-
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
-
-        if (mesh.changing() && 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);
-
-            interface.correct();
-        }
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
-
-        twoPhaseProperties.correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
+
+                if (mesh.changing() && 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);
+
+                    interface.correct();
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H

@@ -1,6 +1,16 @@
-    #include "readTimeControls.H"
+#include "readTimeControls.H"
 
-    bool correctPhi =
+bool correctPhi
-        pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+(
-    bool checkMeshCourantNo =
+    pimple.dict().lookupOrDefault<Switch>("correctPhi", true)
-        pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
+);
+
+bool moveMeshOuterCorrectors
+(
+    pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
+);

applications/solvers/multiphase/interFoam/interFoam.C

@@ -80,14 +80,21 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H

@@ -63,7 +63,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -74,16 +74,23 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        threePhaseProperties.correct();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        #include "alphaEqnsSubCycle.H"
-        interface.correct();
-
-        #define twoPhaseProperties threePhaseProperties
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                threePhaseProperties.correct();
+
+                #include "alphaEqnsSubCycle.H"
+                interface.correct();
+
+                #define twoPhaseProperties threePhaseProperties
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C

@@ -83,14 +83,21 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
+        tmp<surfaceScalarField> tphiAlpha;
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -1,18 +1,4 @@
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("0", dimMass/dimTime, 0)
-);
-
 {
-    #include "alphaControls.H"
-
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
 

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C

@@ -86,44 +86,64 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
-
-        mesh.update();
-
-        if (mesh.changing())
-        {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
-
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
-
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            phi = mesh.Sf() & Uf;
-
-            #include "../interFoam/interDyMFoam/correctPhi.H"
-
-            // Make the flux relative to the mesh motion
-            fvc::makeRelative(phi, U);
-        }
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
-
-        twoPhaseProperties->correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
+
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "../interFoam/interDyMFoam/correctPhi.H"
+
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "alphaControls.H"
+
+            surfaceScalarField rhoPhi
+            (
+                IOobject
+                (
+                    "rhoPhi",
+                    runTime.timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("0", dimMass/dimTime, 0)
+            );
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties->correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -83,14 +83,31 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties->correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            surfaceScalarField rhoPhi
+            (
+                IOobject
+                (
+                    "rhoPhi",
+                    runTime.timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("0", dimMass/dimTime, 0)
+            );
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties->correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -87,14 +87,14 @@ Foam::multiphaseMixture::multiphaseMixture
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
     ),
 
     alphas_

applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H

@@ -54,7 +54,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,

src/finiteVolume/cfdTools/general/include/alphaControls.H

@@ -1,4 +1,12 @@
 const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
 label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+
 label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
-Switch MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
+
+bool MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
+
+bool alphaOuterCorrectors
+(
+    alphaControls.lookupOrDefault<Switch>("alphaOuterCorrectors", false)
+);

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -129,6 +129,9 @@ public:
             //- Helper function to identify when to store the intial residuals
             inline bool storeInitialResiduals() const;
 
+            //- Helper function to identify first PIMPLE (outer) iteration
+            inline bool firstIter() const;
+
             //- Helper function to identify final PIMPLE (outer) iteration
             inline bool finalIter() const;
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,14 +66,20 @@ inline bool Foam::pimpleControl::correct()
 
 inline bool Foam::pimpleControl::storeInitialResiduals() const
 {
-    // start from second PIMPLE iteration
+    // Start from second PIMPLE iteration
     return (corr_ == 2) && (corrPISO_ == 0) && (corrNonOrtho_ == 0);
 }
 
 
+inline bool Foam::pimpleControl::firstIter() const
+{
+    return corr_ == 1;
+}
+
+
 inline bool Foam::pimpleControl::finalIter() const
 {
-    return converged_ || (corr_ == nCorrPIMPLE_);
+    return converged_ || (corr_ == corrPISO_);
 }
 
 

tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict

@@ -47,8 +47,8 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
+maxCo           1;
-maxAlphaCo      0.5;
+maxAlphaCo      1;
 
 maxDeltaT       1;
 

tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution

@@ -19,16 +19,25 @@ solvers
 {
     alpha.water
     {
-        nAlphaCorr      1;
+        nAlphaCorr      2;
-        nAlphaSubCycles 2;
+        nAlphaSubCycles 1;
+        alphaOuterCorrectors yes;
         cAlpha          1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
     pcorr
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-10;
+        tolerance       1e-5;
         relTol          0;
     }
 
@@ -43,14 +52,13 @@ solvers
     p_rghFinal
     {
         $p_rgh;
-        tolerance       1e-07;
         relTol          0;
     }
 
     U
     {
-        solver          PBiCG;
+        solver          smoothSolver;
-        preconditioner  DILU;
+        smoother        symGaussSeidel;
         tolerance       1e-06;
         relTol          0;
     }
@@ -58,10 +66,22 @@ solvers
 
 PIMPLE
 {
-    momentumPredictor no;
+    momentumPredictor   no;
-    nCorrectors     3;
+    nOuterCorrectors    1;
+    nCorrectors         3;
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/polyMesh/boundary

@@ -39,7 +39,7 @@ FoamFile
     bullet
     {
         type            wall;
-        nFaces          37752;
+        nFaces          37743;
         startFace       1133431;
     }
 )