interPhaseChangeFoam: Merged with interPhaseChangeDyMFoam

interPhaseChangeFoam can now run with static or dynamic meshes selected in the
constant/dynamicMeshDict dictionary.

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -106,22 +106,15 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
-                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
-
                 mesh.update();
 
                 if (mesh.changing())
                 {
-
-                    MRF.update();
-
-                    Info<< "Execution time for mesh.update() = "
-                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                        << " s" << endl;
-
                     gh = (g & mesh.C()) - ghRef;
                     ghf = (g & mesh.Cf()) - ghRef;
 
+                    MRF.update();
+
                     if (correctPhi)
                     {
                         // Calculate absolute flux
@@ -143,6 +136,8 @@ int main(int argc, char *argv[])
                 }
             }
 
+            divU = fvc::div(fvc::absolute(phi, U));
+
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 

applications/solvers/multiphase/interFoam/Make/options

@@ -18,6 +18,6 @@ EXE_LIBS = \
     -lincompressibleTurbulenceModels \
     -lfiniteVolume \
     -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
     -lfvOptions \
-    -lmeshTools \
+    -lmeshTools
-    -lsampling

applications/solvers/multiphase/interPhaseChangeFoam/Allwclean

@@ -3,6 +3,5 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 wclean libso phaseChangeTwoPhaseMixtures
 wclean
-wclean interPhaseChangeDyMFoam
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/interPhaseChangeFoam/Allwmake

@@ -6,6 +6,5 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 wmake $targetType phaseChangeTwoPhaseMixtures
 wmake $targetType
-wmake $targetType interPhaseChangeDyMFoam
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/interPhaseChangeFoam/Make/options

@@ -1,14 +1,17 @@
 EXE_INC = \
-    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I. \
+    -I../VoF \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -IphaseChangeTwoPhaseMixtures/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -IphaseChangeTwoPhaseMixtures/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lphaseChangeTwoPhaseMixtures \
@@ -19,6 +22,7 @@ EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lfiniteVolume \
-    -lmeshTools \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
     -lfvOptions \
-    -lsampling
+    -lmeshTools

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,13 +1,16 @@
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
       + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
     );
 
     UEqn.relax();
 
+    fvOptions.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve
@@ -17,10 +20,12 @@
             fvc::reconstruct
             (
                 (
-                    interface.surfaceTensionForce()
+                    mixture.surfaceTensionForce()
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );
+
+        fvOptions.correct(U);
     }

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -2,10 +2,9 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    surfaceScalarField phir("phir", phic*interface.nHatf());
+    surfaceScalarField phir("phir", phic*mixture.nHatf());
 
-    Pair<tmp<volScalarField>> vDotAlphal =
+    Pair<tmp<volScalarField>> vDotAlphal = mixture.vDotAlphal();
-        mixture->vDotAlphal();
     const volScalarField& vDotcAlphal = vDotAlphal[0]();
     const volScalarField& vDotvAlphal = vDotAlphal[1]();
     const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -1,16 +1,14 @@
 {
     // Standard face-flux compression coefficient
-    surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));
+    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
 
     // Add the optional isotropic compression contribution
     if (icAlpha > 0)
     {
         phic *= (1.0 - icAlpha);
-        phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
     }
 
-    volScalarField divU(fvc::div(phi));
-
     if (nAlphaSubCycles > 1)
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();

applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H

@@ -3,7 +3,7 @@ CorrectPhi
     U,
     phi,
     p_rgh,
-    surfaceScalarField("rAUf", fvc::interpolate(rAU)),
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
     divU,
     pimple
 );

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -30,14 +30,18 @@ volVectorField U
 
 
 Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
-autoPtr<phaseChangeTwoPhaseMixture> mixture =
+autoPtr<phaseChangeTwoPhaseMixture> mixturePtr
-    phaseChangeTwoPhaseMixture::New(U, phi);
+(
+    phaseChangeTwoPhaseMixture::New(U, phi)
+);
 
-volScalarField& alpha1(mixture->alpha1());
+phaseChangeTwoPhaseMixture& mixture = mixturePtr();
-volScalarField& alpha2(mixture->alpha2());
 
-const dimensionedScalar& rho1 = mixture->rho1();
+volScalarField& alpha1(mixture.alpha1());
-const dimensionedScalar& rho2 = mixture->rho2();
+volScalarField& alpha2(mixture.alpha2());
+
+const dimensionedScalar& rho1 = mixture.rho1();
+const dimensionedScalar& rho2 = mixture.rho2();
 
 
 // Need to store rho for ddt(rho, U)
@@ -55,13 +59,10 @@ volScalarField rho
 rho.oldTime();
 
 
-// Construct interface from alpha1 distribution
-interfaceProperties interface(alpha1, U, mixture());
-
 // Construct incompressible turbulence model
 autoPtr<incompressible::turbulenceModel> turbulence
 (
-    incompressible::turbulenceModel::New(U, phi, mixture())
+    incompressible::turbulenceModel::New(U, phi, mixture)
 );
 
 

applications/solvers/multiphase/interPhaseChangeFoam/initCorrectPhi.H

@@ -0,0 +1,36 @@
+tmp<volScalarField> rAU;
+
+if (correctPhi)
+{
+    rAU = new volScalarField
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAU", dimTime/dimDensity, 1)
+    );
+
+    volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+
+    #include "correctPhi.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple
+    );
+
+    #include "continuityErrs.H"
+}

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/files

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

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options

@@ -1,30 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I../phaseChangeTwoPhaseMixtures/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lphaseChangeTwoPhaseMixtures \
-    -ltwoPhaseMixture \
-    -linterfaceProperties \
-    -ltwoPhaseProperties \
-    -lincompressibleTransportModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lfiniteVolume \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -lmeshTools \
-    -lfvOptions \
-    -lsampling

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

@@ -1,190 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2018 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
-    interPhaseChangeDyMFoam
-
-Description
-    Solver for 2 incompressible, isothermal immiscible fluids with phase-change
-    (e.g. cavitation).  Uses a VOF (volume of fluid) phase-fraction based
-    interface capturing approach, with optional mesh motion and mesh topology
-    changes including adaptive re-meshing.
-
-    The momentum and other fluid properties are of the "mixture" and a
-    single momentum equation is solved.
-
-    The set of phase-change models provided are designed to simulate cavitation
-    but other mechanisms of phase-change are supported within this solver
-    framework.
-
-    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "CMULES.H"
-#include "subCycle.H"
-#include "interfaceProperties.H"
-#include "phaseChangeTwoPhaseMixture.H"
-#include "turbulentTransportModel.H"
-#include "pimpleControl.H"
-#include "fvOptions.H"
-#include "CorrectPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "initContinuityErrs.H"
-    #include "createFields.H"
-
-    volScalarField rAU
-    (
-        IOobject
-        (
-            "rAU",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
-    );
-
-    #include "createUf.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    turbulence->validate();
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readDyMControls.H"
-
-        // Store divU from the previous mesh so that it can be mapped
-        // and used in correctPhi to ensure the corrected phi has the
-        // same divergence
-        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
-
-        #include "CourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // --- 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()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-                }
-
-                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);
-                }
-
-                if (mesh.changing() && checkMeshCourantNo)
-                {
-                    #include "meshCourantNo.H"
-                }
-            }
-
-            #include "alphaControls.H"
-
-            surfaceScalarField rhoPhi
-            (
-                IOobject
-                (
-                    "rhoPhi",
-                    runTime.timeName(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar("0", dimMass/dimTime, 0)
-            );
-
-            mixture->correct();
-
-            #include "alphaEqnSubCycle.H"
-            interface.correct();
-
-            #include "UEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H

@@ -1,80 +0,0 @@
-{
-    rAU = 1.0/UEqn.A();
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(HbyA)
-      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
-    );
-
-    if (p_rgh.needReference())
-    {
-        fvc::makeRelative(phiHbyA, U);
-        adjustPhi(phiHbyA, U, p_rgh);
-        fvc::makeAbsolute(phiHbyA, U);
-    }
-
-    surfaceScalarField phig
-    (
-        (
-            interface.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
-    );
-
-    phiHbyA += phig;
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
-
-    Pair<tmp<volScalarField>> vDotP = mixture->vDotP();
-    const volScalarField& vDotcP = vDotP[0]();
-    const volScalarField& vDotvP = vDotP[1]();
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqn
-        (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-          - (vDotvP - vDotcP)*(mixture->pSat() - rho*gh)
-          + fvm::Sp(vDotvP - vDotcP, p_rgh)
-        );
-
-        p_rghEqn.setReference(pRefCell, pRefValue);
-
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + p_rghEqn.flux();
-
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-        }
-    }
-
-    {
-        Uf = fvc::interpolate(U);
-        surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += n*(phi/mesh.magSf() - (n & Uf));
-    }
-
-    // Make the fluxes relative to the mesh motion
-    fvc::makeRelative(phi, U);
-
-    p == p_rgh + rho*gh;
-
-    if (p_rgh.needReference())
-    {
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
-        );
-        p_rgh = p - rho*gh;
-    }
-}

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -27,7 +27,8 @@ Application
 Description
     Solver for 2 incompressible, isothermal immiscible fluids with phase-change
     (e.g. cavitation).  Uses a VOF (volume of fluid) phase-fraction based
-    interface capturing approach.
+    interface capturing approach, with optional mesh motion and mesh topology
+    changes including adaptive re-meshing.
 
     The momentum and other fluid properties are of the "mixture" and a
     single momentum equation is solved.
@@ -41,6 +42,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
 #include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
@@ -48,6 +50,7 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
+#include "CorrectPhi.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,23 +60,33 @@ int main(int argc, char *argv[])
 
     #include "setRootCaseLists.H"
     #include "createTime.H"
-    #include "createMesh.H"
+    #include "createDynamicFvMesh.H"
-    #include "createControl.H"
+    #include "createDyMControls.H"
+    #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "createTimeControls.H"
+    #include "initCorrectPhi.H"
-    #include "CourantNo.H"
+    #include "createUfIfPresent.H"
-    #include "setInitialDeltaT.H"
 
     turbulence->validate();
 
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "readDyMControls.H"
+
+        // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
@@ -83,7 +96,37 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            #include "alphaControls.H"
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    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);
+
+                        mixture.correct();
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
+                }
+            }
+
+            divU = fvc::div(fvc::absolute(phi, U));
 
             surfaceScalarField rhoPhi
             (
@@ -97,10 +140,10 @@ int main(int argc, char *argv[])
                 dimensionedScalar("0", dimMass/dimTime, 0)
             );
 
-            mixture->correct();
+            #include "alphaControls.H"
-
             #include "alphaEqnSubCycle.H"
-            interface.correct();
+
+            mixture.correct();
 
             #include "UEqn.H"
 

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -1,19 +1,33 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn.A());
+    if (correctPhi)
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+    {
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+    volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         fvc::flux(HbyA)
-      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+      + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)
     );
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
         adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
 
     surfaceScalarField phig
     (
         (
-            interface.surfaceTensionForce()
+            mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf()
     );
@@ -23,7 +37,7 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, U, phiHbyA, rAUf);
 
-    Pair<tmp<volScalarField>> vDotP = mixture->vDotP();
+    Pair<tmp<volScalarField>> vDotP = mixture.vDotP();
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();
 
@@ -32,7 +46,7 @@
         fvScalarMatrix p_rghEqn
         (
             fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-          - (vDotvP - vDotcP)*(mixture->pSat() - rho*gh)
+          - (vDotvP - vDotcP)*(mixture.pSat() - rho*gh)
           + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -44,12 +58,18 @@
         {
             phi = phiHbyA + p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
     }
 
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
     p == p_rgh + rho*gh;
 
     if (p_rgh.needReference())
@@ -62,4 +82,9 @@
         );
         p_rgh = p - rho*gh;
     }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
 }

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C

@@ -66,8 +66,8 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::Kunz
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
 
     return Pair<tmp<volScalarField>>
     (
@@ -81,8 +81,8 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
 
     return Pair<tmp<volScalarField>>
     (

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options

@@ -1,11 +1,11 @@
 EXE_INC = \
-    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
-    -ltwoPhaseMixture \
+    -limmiscibleIncompressibleTwoPhaseMixture \
-    -ltwoPhaseProperties \
-    -lincompressibleTransportModels \
     -lfiniteVolume

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C

@@ -66,7 +66,7 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::Merkle
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotAlphal() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
 
     return Pair<tmp<volScalarField>>
     (
@@ -78,8 +78,8 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotAlphal() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
 
     return Pair<tmp<volScalarField>>
     (

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

@@ -96,7 +96,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff
     const volScalarField& p
 ) const
 {
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
     volScalarField rho
     (
         limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
@@ -111,10 +111,10 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
     volScalarField pCoeff(this->pCoeff(p));
 
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
 
     return Pair<tmp<volScalarField>>
     (
@@ -128,10 +128,10 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
 {
-    const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
+    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
     volScalarField pCoeff(this->pCoeff(p));
 
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
     volScalarField apCoeff(limitedAlpha1*pCoeff);
 
     return Pair<tmp<volScalarField>>

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C

@@ -42,7 +42,7 @@ Foam::phaseChangeTwoPhaseMixture::phaseChangeTwoPhaseMixture
     const surfaceScalarField& phi
 )
 :
-    incompressibleTwoPhaseMixture(U, phi),
+    immiscibleIncompressibleTwoPhaseMixture(U, phi),
     phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs")),
     pSat_("pSat", dimPressure, lookup("pSat"))
 {}
@@ -53,7 +53,7 @@ Foam::phaseChangeTwoPhaseMixture::phaseChangeTwoPhaseMixture
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
 {
-    volScalarField alphalCoeff(1.0/rho1() - alpha1_*(1.0/rho1() - 1.0/rho2()));
+    volScalarField alphalCoeff(1.0/rho1() - alpha1()*(1.0/rho1() - 1.0/rho2()));
     Pair<tmp<volScalarField>> mDotAlphal = this->mDotAlphal();
 
     return Pair<tmp<volScalarField>>
@@ -75,7 +75,7 @@ Foam::phaseChangeTwoPhaseMixture::vDotP() const
 
 bool Foam::phaseChangeTwoPhaseMixture::read()
 {
-    if (incompressibleTwoPhaseMixture::read())
+    if (immiscibleIncompressibleTwoPhaseMixture::read())
     {
         phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
         lookup("pSat") >> pSat_;

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H

@@ -35,7 +35,7 @@ SourceFiles
 #ifndef phaseChangeTwoPhaseMixture_H
 #define phaseChangeTwoPhaseMixture_H
 
-#include "incompressibleTwoPhaseMixture.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
 #include "typeInfo.H"
 #include "runTimeSelectionTables.H"
 #include "volFields.H"
@@ -54,7 +54,7 @@ namespace Foam
 
 class phaseChangeTwoPhaseMixture
 :
-    public incompressibleTwoPhaseMixture
+    public immiscibleIncompressibleTwoPhaseMixture
 {
 
 protected:

bin/interPhaseChangeDyMFoam

@@ -0,0 +1 @@
+mergedDyM

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict

@@ -32,7 +32,7 @@ writeInterval           0.001;
 
 purgeWrite              0;
 
-writeFormat             ascii;
+writeFormat             binary;
 
 writePrecision          6;
 
@@ -45,6 +45,6 @@ runTimeModifiable       yes;
 adjustTimeStep          on;
 
 maxCo                   5;
-
+maxAlphaCo              2;
 
 // ************************************************************************* //

tutorials/multiphase/interPhaseChangeFoam/propeller/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     interPhaseChangeDyMFoam;
+application     interPhaseChangeFoam;
 
 startFrom       startTime;