INT: Initial update of isoAdvector and interIsoFoam to work with AMR.

applications/solvers/multiphase/interIsoFoam/Make/options

@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
@@ -15,6 +16,7 @@ EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lfiniteVolume \
+    -ldynamicFvMesh \
     -lfvOptions \
     -lmeshTools \
     -lsampling \

applications/solvers/multiphase/interIsoFoam/alphaControls.H

@@ -1,3 +1,29 @@
 const dictionary& alphaControls = mesh.solverDict(alpha1.name());
 
+//label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+
 label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
+/*
+bool MULESCorr(alphaControls.lookupOrDefault("MULESCorr", false));
+
+// Apply the compression correction from the previous iteration
+// Improves efficiency for steady-simulations but can only be applied
+// once the alpha field is reasonably steady, i.e. fully developed
+bool alphaApplyPrevCorr
+(
+    alphaControls.lookupOrDefault("alphaApplyPrevCorr", false)
+);
+
+// Isotropic compression coefficient
+scalar icAlpha
+(
+    alphaControls.lookupOrDefault<scalar>("icAlpha", 0)
+);
+
+// Shear compression coefficient
+scalar scAlpha
+(
+    alphaControls.lookupOrDefault<scalar>("scAlpha", 0)
+);
+*/

applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H

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

applications/solvers/multiphase/interIsoFoam/alphaEqn.H

@@ -1,35 +1,50 @@
-// If there are more than one outer corrector, we use a mixture of old and
-// new U and phi for propagating alpha1 in all but the first outer iteration
-if (!pimple.firstIter())
-{
-    // We are recalculating alpha1 from the its old time value
-    alpha1 = alpha1.oldTime();
-    // Temporarily storing new U and phi values in prevIter storage
-    U.storePrevIter();
-    phi.storePrevIter();
+    {
+        // Note for AMR implementation:
+        // At this point we have just entered the new time step,
+        // the mesh has been refined and the alpha, phi and U contain
+        // the field values at the beginning of this time step mapped
+        // to the new mesh.
 
-    // Overwriting new U and phi values with mixture of old and new values
-    phi = 0.5*(phi + phi.oldTime());
-    U = 0.5*(U + U.oldTime());
-}
+        // The above assumes that we are in firstIter() of the outer
+        // corrector loop. If we are in any subsequent iter of this loop
+        // the alpha1, U and phi will be overwritten with the new time step
+        // values but still on the same mesh.
 
-// Update alpha1
-advector.advect();
 
-// Update rhoPhi
-rhoPhi = advector.getRhoPhi(rho1, rho2);
+        if (pimple.firstIter())
+        {
+            // Note: This assumes moveMeshOuterCorrectors is false
+            alpha10 = alpha1;
+            U0 = U;
+            phi0 = phi;
+            advector.advect();
 
-alpha2 = 1.0 - alpha1;
+            #include "rhofs.H"
+            rhoPhi = advector.getRhoPhi(rho1f, rho2f);
+        }
+        else
+        {
+            alpha1 = alpha10;
+            // Temporarily setting U and phi to average of old and new value
+            // Note: Illegal additions if mesh is topoChanging
+            // (e.g. if moveMeshOuterCorrectors and AMR)
+            U = 0.5*U0 + 0.5*U;
+            phi = 0.5*phi0 + 0.5*phi;
+            isoAdvection advector(alpha1, phi, U);
+            advector.advect();
+            #include "rhofs.H"
+            rhoPhi = advector.getRhoPhi(rho1f, rho2f);
+            // Resetting U and phi to the new value
+            U = 2.0*U - U0;
+            phi = 2.0*phi - phi0;
+        }
 
-if (!pimple.firstIter())
-{
-    // Restoring new U and phi values temporarily saved in prevIter() above
-    U = U.prevIter();
-    phi = phi.prevIter();
-}
+        alpha2 = 1.0 - alpha1;
+        mixture.correct();
 
-Info<< "Phase-1 volume fraction = "
-    << alpha1.weightedAverage(mesh.Vsc()).value()
-    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
-    << endl;
+    }
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+        << endl;

applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H

@@ -15,6 +15,12 @@ if (nAlphaSubCycles > 1)
 
     tmp<volScalarField> trSubDeltaT;
 
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
     for
     (
         subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);

applications/solvers/multiphase/interIsoFoam/alphaSuSp.H

@@ -0,0 +1,3 @@
+zeroField Su;
+zeroField Sp;
+zeroField divU;

applications/solvers/multiphase/interIsoFoam/correctPhi.H

@@ -3,7 +3,7 @@ CorrectPhi
     U,
     phi,
     p_rgh,
-    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
     geometricZeroField(),
     pimple
 );

applications/solvers/multiphase/interIsoFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (
@@ -119,5 +121,67 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
+/*
+// MULES compressed flux is registered in case scalarTransport FO needs it.
+surfaceScalarField alphaPhiUn
+(
+    IOobject
+    (
+        "alphaPhiUn",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(phi.dimensions(), Zero)
+);
+*/
+
 #include "createMRF.H"
-#include "createIsoAdvection.H"
+#include "createFvOptions.H"
+
+//Auxiliary fields when nOuterCorrectors > 1
+surfaceScalarField phi0
+(
+    IOobject
+    (
+        "phi0",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    phi
+);
+
+
+volVectorField U0
+(
+    IOobject
+    (
+        "U0",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    U
+);
+
+
+volScalarField alpha10
+(
+    IOobject
+    (
+        "alpha10",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    alpha1
+);
+
+
+isoAdvection advector(alpha1, phi, U);

applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H

@@ -1,2 +0,0 @@
-// Defining isoAdvection
-isoAdvection advector(alpha1, phi, U);

applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H

@@ -0,0 +1,34 @@
+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)
+    );
+
+    #include "correctPhi.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple
+    );
+
+    #include "continuityErrs.H"
+}

applications/solvers/multiphase/interIsoFoam/interIsoFoam.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,36 +32,35 @@ Group
 
 Description
     Solver derived from interFoam for 2 incompressible, isothermal immiscible
-    fluids using the iso-advector phase-fraction based interface capturing
-    approach.
-
-    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.
-
-    For a two-fluid approach see twoPhaseEulerFoam.
+    fluids using the isoAdvector phase-fraction based interface capturing
+    approach, with optional mesh motion and mesh topology changes including
+    adaptive re-meshing.
 
     Reference:
-        \verbatim
-            Roenby, J., Bredmose, H. and Jasak, H. (2016).
-            A computational method for sharp interface advection
-            Royal Society Open Science, 3
-            doi 10.1098/rsos.160405
-        \endverbatim
+    \verbatim
+        Roenby, J., Bredmose, H. and Jasak, H. (2016).
+        A computational method for sharp interface advection
+        Royal Society Open Science, 3
+        doi 10.1098/rsos.160405
+    \endverbatim
 
-    isoAdvector code supplied by Johan Roenby, DHI (2016)
+    isoAdvector code supplied by Johan Roenby, STROMNING (2018)
 
 \*---------------------------------------------------------------------------*/
 
-#include "isoAdvection.H"
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "isoAdvection.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "immiscibleIncompressibleTwoPhaseMixture.H"
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,31 +71,39 @@ int main(int argc, char *argv[])
     #include "addCheckCaseOptions.H"
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
+    #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
-    #include "correctPhi.H"
+//    #include "createAlphaFluxes.H"
+    #include "initCorrectPhi.H"
+    #include "createUfIfPresent.H"
 
     turbulence->validate();
 
-    #include "readTimeControls.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
+    if (!LTS)
+    {
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "readDyMControls.H"
 
-        #include "CourantNo.H"
-        #include "alphaCourantNo.H"
-        #include "setDeltaT.H"
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
 
@@ -104,6 +112,45 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    // Do not apply previous time-step mesh compression flux
+                    // if the mesh topology changed
+                    if (mesh.topoChanging())
+                    {
+//                        talphaPhi1Corr0.clear();
+                    }
+
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    MRF.update();
+
+                    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"
+                    }
+                }
+            }
+
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 

applications/solvers/multiphase/interIsoFoam/pEqn.H

@@ -1,17 +1,29 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    if (correctPhi)
+    {
+        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)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
     );
     MRF.makeRelative(phiHbyA);
-    adjustPhi(phiHbyA, U, p_rgh);
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
 
     surfaceScalarField phig
     (
@@ -19,7 +31,6 @@
             mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf()
-//          - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf
     );
 
     phiHbyA += phig;
@@ -44,7 +55,7 @@
 
             p_rgh.relax();
 
-            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);
         }
@@ -52,6 +63,12 @@
 
     #include "continuityErrs.H"
 
+    // 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())
@@ -64,4 +81,9 @@
         );
         p_rgh = p - rho*gh;
     }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
 }

applications/solvers/multiphase/interIsoFoam/rhofs.H

@@ -0,0 +1,2 @@
+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);

applications/solvers/multiphase/interIsoFoam/setDeltaT.H

@@ -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-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/multiphase/interIsoFoam/setRDeltaT.H

@@ -0,0 +1,136 @@
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    scalar maxCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    scalar alphaSpreadDiff
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
+    );
+
+    scalar alphaSpreadMax
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    scalar alphaSpreadMin
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    label nAlphaSweepIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.ref() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(rhoPhi))()()
+       /((2*maxCo)*mesh.V()*rho())
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        volScalarField alpha1Bar(fvc::average(alpha1));
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            pos0(alpha1Bar() - alphaSpreadMin)
+           *pos0(alphaSpreadMax - alpha1Bar())
+           *fvc::surfaceSum(mag(phi))()()
+           /((2*maxAlphaCo)*mesh.V())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread
+        (
+            rDeltaT,
+            alpha1,
+            nAlphaSpreadIter,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.primitiveField())
+            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+}