ENH: Adding functionality to scalarTransport FO and residence time tutorials for VOF

and single phase cases. Registration of the compressed flux in interFoam as it is
needed for the FO if used.

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -126,35 +126,35 @@
     {
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
+        alphaPhiUn =
             (
-               -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha1,
-                alpharScheme
-            )
-        );
+                fvc::flux
+                (
+                    phi,
+                    alpha1,
+                    alphaScheme
+                )
+            + fvc::flux
+                (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+            );
 
         // Calculate the Crank-Nicolson off-centred alpha flux
         if (ocCoeff > 0)
         {
-            talphaPhiUn =
-                cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
+            alphaPhiUn =
+                cnCoeff*alphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
         }
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhiCorr(alphaPhiUn - alphaPhi);
             volScalarField alpha10("alpha10", alpha1);
 
-            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);
+            MULES::correct(alpha1, alphaPhiUn, talphaPhiCorr.ref(), 1, 0);
 
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
@@ -169,7 +169,7 @@
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi = alphaPhiUn;
 
             MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
         }

applications/solvers/multiphase/interFoam/createFields.H

@@ -135,6 +135,21 @@ surfaceScalarField alphaPhi
     phi*fvc::interpolate(alpha1)
 );
 
+// 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("zero", phi.dimensions(), 0.0)
+);
+
 // MULES Correction
 tmp<surfaceScalarField> talphaPhiCorr0;
 

src/functionObjects/solvers/scalarTransport/scalarTransport.C

@@ -29,6 +29,7 @@ License
 #include "fvmDiv.H"
 #include "fvmLaplacian.H"
 #include "fvmSup.H"
+#include "CMULES.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
 #include "addToRunTimeSelectionTable.H"
@@ -85,7 +86,8 @@ Foam::volScalarField& Foam::functionObjects::scalarTransport::transportedField()
 Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
 (
     const volScalarField& s,
-    const surfaceScalarField& phi
+    const surfaceScalarField& phi,
+    const volScalarField& alpha
 ) const
 {
     typedef incompressible::turbulenceModel icoModel;
@@ -93,9 +95,11 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
 
     word Dname("D" + s.name());
 
+    volScalarField phaseMask(pos(alpha - 0.99));
+
     if (constantD_)
     {
-        return tmp<volScalarField>
+        tmp<volScalarField> tD
         (
             new volScalarField
             (
@@ -111,6 +115,18 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
                 dimensionedScalar(Dname, phi.dimensions()/dimLength, D_)
             )
         );
+
+        return phaseMask*tD;
+    }
+    else if (nutName_ != "none")
+    {
+        const volScalarField& nutMean =
+            mesh_.lookupObject<volScalarField>(nutName_);
+
+        return tmp<volScalarField>
+        (
+            new volScalarField(Dname, phaseMask*nutMean)
+        );
     }
     else if (foundObject<icoModel>(turbulenceModel::propertiesName))
     {
@@ -119,7 +135,10 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             turbulenceModel::propertiesName
         );
 
-        return model.nuEff();
+        return tmp<volScalarField>
+        (
+             new volScalarField(Dname, phaseMask*model.nuEff())
+        );
     }
     else if (foundObject<cmpModel>(turbulenceModel::propertiesName))
     {
@@ -128,7 +147,10 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             turbulenceModel::propertiesName
         );
 
-        return model.muEff();
+        return tmp<volScalarField>
+        (
+             new volScalarField(Dname, phaseMask*model.muEff())
+        );
     }
     else
     {
@@ -163,14 +185,22 @@ Foam::functionObjects::scalarTransport::scalarTransport
 :
     fvMeshFunctionObject(name, runTime, dict),
     fieldName_(dict.lookupOrDefault<word>("field", "s")),
-    phiName_("phi"),
-    rhoName_("rho"),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    UPhiName_(dict.lookupOrDefault<word>("UPhi", "none")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    nutName_(dict.lookupOrDefault<word>("nut", "none")),
+    phaseName_(dict.lookupOrDefault<word>("phase", "none")),
+    phasePhiCompressedName_
+    (
+        dict.lookupOrDefault<word>("phasePhiCompressed", "alphaPhiUn")
+    ),
     D_(0),
     constantD_(false),
     nCorr_(0),
     resetOnStartUp_(false),
     schemesField_("unknown-schemesField"),
-    fvOptions_(mesh_)
+    fvOptions_(mesh_),
+    bounded01_(dict.lookupOrDefault<bool>("bounded01", true))
 {
     read(dict);
 
@@ -199,6 +229,11 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
 
     dict.readIfPresent("phi", phiName_);
     dict.readIfPresent("rho", rhoName_);
+    dict.readIfPresent("UPhi", UPhiName_);
+    dict.readIfPresent("nut", nutName_);
+    dict.readIfPresent("phase", phaseName_);
+    dict.readIfPresent("bounded01", bounded01_);
+
     schemesField_ = dict.lookupOrDefault("schemesField", fieldName_);
 
     constantD_ = false;
@@ -223,12 +258,93 @@ bool Foam::functionObjects::scalarTransport::execute()
 {
     Log << type() << " write:" << endl;
 
-    const surfaceScalarField& phi = lookupObject<surfaceScalarField>(phiName_);
+    tmp<surfaceScalarField> tPhi
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "phi",
+                mesh_.time().timeName(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("tPhi", dimMass/dimTime, 0.0)
+        )
+    );
+    surfaceScalarField& phi = tPhi.ref();
+
+    const dimensionSet dim
+    (
+        mesh_.lookupObject<surfaceScalarField>(phiName_).dimensions()
+    );
+
+    bool compressible = true;
+    if (dim == dimVolume/dimTime)
+    {
+        compressible = false;
+        phi.dimensions().reset(dimVolume/dimTime);
+    }
+
+    //Obtain phi from phiName or constructed from UPhiName
+    if (phiName_ != "none")
+    {
+        phi = const_cast<surfaceScalarField&>
+        (
+            mesh_.lookupObject<surfaceScalarField>(phiName_)
+        );
+    }
+    else if(UPhiName_ != "none")
+    {
+        const volVectorField& Uphi =
+            mesh_.lookupObject<volVectorField>(UPhiName_);
+
+        if (!compressible)
+        {
+            phi = fvc::interpolate(Uphi) & mesh_.Sf();
+        }
+        else
+        {
+            const volScalarField& rho =
+                mesh_.lookupObject<volScalarField>(rhoName_);
+
+            phi = fvc::interpolate(rho*Uphi) & mesh_.Sf();
+        }
+    }
+
+    tmp<volScalarField> tPhaseMask
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "tPhaseMask",
+                mesh_.time().timeName(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("tPhaseMask", dimless, 1.0)
+        )
+    );
+    volScalarField& phaseMask = tPhaseMask.ref();
+
+    // Set phaseMask if s is transported in a phase
+    if (phaseName_ != "none")
+    {
+        const volScalarField& alpha =
+            mesh_.lookupObject<volScalarField>(phaseName_);
+
+        phaseMask = alpha;
+    }
 
     volScalarField& s = transportedField();
 
     // Calculate the diffusivity
-    volScalarField D(this->D(s, phi));
+    volScalarField D(this->D(s, phi, phaseMask));
 
     word divScheme("div(phi," + schemesField_ + ")");
     word laplacianScheme("laplacian(" + D.name() + "," + schemesField_ + ")");
@@ -240,12 +356,84 @@ bool Foam::functionObjects::scalarTransport::execute()
         relaxCoeff = mesh_.equationRelaxationFactor(schemesField_);
     }
 
-    if (phi.dimensions() == dimMass/dimTime)
+    // two phase scalar transport
+    if (phaseName_ != "none")
+    {
+        const volScalarField& alpha =
+            mesh_.lookupObject<volScalarField>(phaseName_);
+
+        const surfaceScalarField& limitedPhiAlpa =
+            mesh_.lookupObject<surfaceScalarField>(phasePhiCompressedName_);
+
+/*
+        surfaceScalarField phic(2.0*mag(phi/mesh_.magSf()));
+
+        const volVectorField gradAlpha(fvc::grad(alpha, "nHat"));
+
+        surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
+
+        dimensionedScalar deltaN
+        (
+            "deltaN", 1e-8/pow(average(mesh_.V()), 1.0/3.0)
+        );
+
+        surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN));
+
+        surfaceScalarField nHat(nHatfv & mesh_.Sf());
+
+        surfaceScalarField phir(phic*nHat);
+
+        surfaceScalarField limitedPhiAlpa
+        (
+            fvc::flux
+            (
+                phi,
+                alpha,
+                "div(phi,s)"
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, (1-alpha), "div(phirb,s)"),
+                alpha,
+                "div(phirb,s)"
+            )
+        );
+*/
+        // Reset D dimensions consistent with limitedPhiAlpa
+        D.dimensions().reset(limitedPhiAlpa.dimensions()/dimLength);
+
+        tmp<surfaceScalarField> tTPhiUD;
+        // Solve
+        for (label i = 0; i <= nCorr_; i++)
+        {
+            fvScalarMatrix sEqn
+            (
+                fvm::ddt(s)
+              + fvm::div(limitedPhiAlpa, s, divScheme)
+              - fvm::laplacian(D, s, laplacianScheme)
+              ==
+                alpha*fvOptions_(s)
+            );
+
+            sEqn.relax(relaxCoeff);
+            fvOptions_.constrain(sEqn);
+            sEqn.solve(mesh_.solverDict(schemesField_));
+
+            tTPhiUD = sEqn.flux();
+        }
+
+        if (bounded01_)
+        {
+            MULES::explicitSolve(s, phi, tTPhiUD.ref(), 1, 0);
+        }
+    }
+    else if (compressible)
     {
         const volScalarField& rho = lookupObject<volScalarField>(rhoName_);
 
         for (label i = 0; i <= nCorr_; i++)
         {
+
             fvScalarMatrix sEqn
             (
                 fvm::ddt(rho, s)
@@ -262,7 +450,7 @@ bool Foam::functionObjects::scalarTransport::execute()
             sEqn.solve(mesh_.solverDict(schemesField_));
         }
     }
-    else if (phi.dimensions() == dimVolume/dimTime)
+    else if (!compressible)
     {
         for (label i = 0; i <= nCorr_; i++)
         {

src/functionObjects/solvers/scalarTransport/scalarTransport.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,8 +33,12 @@ Description
     - To specify the field name set the 'field' entry
     - To employ the same numerical schemes as another field set
       the 'schemesField' entry,
-    - The diffusivity can be set manually using the 'D' entry, or retrieved
-      from the turbulence model (if applicable).
+    - The diffusivity can be set manually using the 'D' entry, retrieved
+      from the turbulence model or specified nut
+    - To specify a different flux derived from U enter UPhi velocity field name
+      (the phi used will be calculated based on this UPhi)
+    - To specify a transport quantity within a phase enter phase.
+    - bounded01 bounds the transported scalar within 0 and 1.
 
 Usage
     Example of function object specification to solve a scalar transport
@@ -47,6 +51,11 @@ Usage
             type            scalarTransport;
             libs            ("libutilityFunctionObjects.so");
 
+            resetOnStartUp  no;
+            region          cabin;
+            field           H2O;
+
+
             fvOptions
             {
                 ...
@@ -55,17 +64,67 @@ Usage
     }
     \endverbatim
 
+    Example of function object specification to solve a residency time
+    in a two phase flow:
+    equation:
+    \verbatim
+    functions
+    {
+        sTransport
+        {
+            type            scalarTransport;
+            libs            ("libsolverFunctionObjects.so");
+
+            enabled         true;
+            writeControl    outputTime;
+            writeInterval   1;
+
+            field           s;
+            bounded01       false;
+            phase           alpha.water;
+
+            write           true;
+
+            fvOptions
+            {
+                unitySource
+                {
+                    type            scalarSemiImplicitSource;
+                    enabled         true;
+
+                    scalarSemiImplicitSourceCoeffs
+                    {
+                        selectionMode   all;
+                        volumeMode      specific;
+                        injectionRateSuSp
+                        {
+                            s           (1 0);
+                        }
+                    }
+                }
+            }
+
+            resetOnStartUp  false;
+        }
+    }
+    \endverbatim
+
     Where the entries comprise:
     \table
         Property     | Description             | Required    | Default value
         type         | Type name: scalarTransport | yes      |
         phi          | Name of flux field      | no          | phi
         rho          | Name of density field   | no          | rho
-        D            | Diffision coefficient   | no          | auto generated
+        UPhi         | Name of U to generate phi  | no          | none
+        phase        | Name of the phase name  | no          | none
+        nut          | Name of the turbulence viscocity | no | none
+        D            | Diffusion coefficient   | no          | auto generated
         nCorr        | Number of correctors    | no          | 0
         resetOnStartUp | Reset scalar to zero on start-up | no | no
         schemesField | Name of field to specify schemes | no | fieldName
         fvOptions    | List of scalar sources  | no          |
+        bounded01    | Bounds scalar betwee 0-1 for multiphase | no |true
+        phasePhiCompressed |Compressed flux for VOF | no | alphaPhiUn
     \endtable
 
 See also
@@ -106,9 +165,22 @@ class scalarTransport
         //- Name of flux field (optional)
         word phiName_;
 
+        //- Name of velocity field from which the flux is obtained if phiName is
+        //  not given (optional)
+        word UPhiName_;
+
         //- Name of density field (optional)
         word rhoName_;
 
+        //- Name of turbulent viscosity field (optional)
+        word nutName_;
+
+        //- Name of phase field
+        word phaseName_;
+
+        //- Name of phase field compressed flux
+        word phasePhiCompressedName_;
+
         //- Diffusion coefficient (optional)
         scalar D_;
 
@@ -127,6 +199,9 @@ class scalarTransport
         //- Run-time selectable finite volume options, e.g. sources, constraints
         fv::optionList fvOptions_;
 
+        //- Bound scalar between 0-1 using MULES for multiphase case
+        bool bounded01_;
+
 
     // Private Member Functions
 
@@ -137,7 +212,8 @@ class scalarTransport
         tmp<volScalarField> D
         (
             const volScalarField& s,
-            const surfaceScalarField& phi
+            const surfaceScalarField& phi,
+            const volScalarField& alpha
         ) const;
 
         //- Disallow default bitwise copy construct

tutorials/incompressible/pimpleFoam/RAS/TJunction/0/s

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      s;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet1
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+
+    outlet2
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/RAS/TJunction/system/controlDict

@@ -79,6 +79,41 @@ functions
             (0.21 0 0.01)         // at central block
         );
     }
+
+    sTransport
+    {
+        type            scalarTransport;
+        libs            ("libsolverFunctionObjects.so");
+
+        enabled         true;
+        writeControl    outputTime;
+        writeInterval   1;
+
+        field           s;
+
+        write           true;
+
+        fvOptions
+        {
+            unitySource
+            {
+                type            scalarSemiImplicitSource;
+                enabled         true;
+
+                scalarSemiImplicitSourceCoeffs
+                {
+                    selectionMode   all;