Merge remote-tracking branch 'origin/develop'

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "cfmesh"]
+	path = modules/cfmesh
+	url = https://develop.openfoam.com/Community/integration-cfmesh.git
+[submodule "avalanche"]
+	path = modules/avalanche
+	url = https://develop.openfoam.com/Community/avalanche.git

Allwmake

@@ -43,6 +43,15 @@ echo "Compile OpenFOAM applications"
 echo
 applications/Allwmake $targetType $*
 
+# Additional components/modules
+if [ -d "$WM_PROJECT_DIR/modules" ]
+then
+    echo "========================================"
+    echo "Compile OpenFOAM modules"
+    echo
+    (cd $WM_PROJECT_DIR/modules 2>/dev/null && wmake -all)
+fi
+
 # Some summary information
 echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"

BuildIssues.txt

@@ -1,8 +1,23 @@
-OpenFOAM-1706
+OpenFOAM-1712
 ==================
 Known Build Issues
 ==================
 
+---------------------
+Intel MPI (Gcc/Clang)
+---------------------
+
+  Either I_MPI_ROOT or MPI_ROOT can be used to specify the Intel-MPI
+  installation directory path.
+
+  The ThirdParty build of ptscotch uses `mpiicc` for Intel-MPI
+  instead of the usual `mpicc`.
+  When gcc or clang are used, it is highly likely that the
+  I_MPI_CC environment variable also needs to be set accordingly.
+
+  See `mpiicc -help` for more information about environment variables.
+
+
 --------------
 Intel Compiler
 --------------
@@ -60,6 +75,33 @@ If your system compiler is too old to build the minimum required gcc or
 clang/llvm, it is just simply too old.
 
 
+---------------------------------
+ThirdParty clang without gmp/mpfr
+---------------------------------
+
+If using ThirdParty clang without gmp/mpfr, the ThirdParty makeCGAL
+script will need to be run manually and specify that there is no
+gmp/mpfr. Eg,
+
+    cd $WM_THIRD_PARTY_DIR
+    ./makeCGAL gmp-none mpfr-none
+
+Subequent compilation with Allwmake will now run largely without any
+problems, except that the components linking against CGAL
+(foamyMesh and surfaceBooleanFeatures) will also try to link against
+a nonexistent mpfr library. As a workaround, the link-dependency can
+be removed in wmake/rules/General/CGAL :
+
+    CGAL_LIBS = \
+        -L$(BOOST_ARCH_PATH)/lib \
+        -L$(BOOST_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) \
+        -L$(CGAL_ARCH_PATH)/lib \
+        -L$(CGAL_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) \
+        -lCGAL
+
+This is a temporary inconvenience until a more robust solution is found.
+
+
 -------------------------
 Building with spack
 -------------------------

README.md

@@ -1,7 +1,8 @@
 # About OpenFOAM
 OpenFOAM is a free, open source CFD software [released and developed primarily by OpenCFD Ltd](http://www.openfoam.com) since 2004. It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.  [More...](http://www.openfoam.com/documentation)
 
-OpenFOAM+ is professionally released every six months to include customer sponsored developments and contributions from the community, including the OpenFOAM Foundation. Releases designated OpenFOAM+ contain several man years of client-sponsored developments of which much has been transferred to, but not released in the OpenFOAM Foundation branch.
+
+OpenFOAM is professionally released every six months to include customer sponsored developments and contributions from the community - individual and group contributors, fork re-integrations including from FOAM-extend and OpenFOAM Foundation Ltd - in this Official Release sanctioned by the OpenFOAM Worldwide Trademark Owner aiming towards one OpenFOAM.
 
 
 # Copyright
@@ -9,7 +10,7 @@ OpenFOAM is free software: you can redistribute it and/or modify it under the te
 
 
 # OpenFOAM Trademark
-OpenCFD Ltd grants use of the OpenFOAM trademark by Third Parties on a licence basis. ESI Group and the OpenFOAM Foundation Ltd are currently permitted to use the Name and agreed Domain Name. For information on trademark use, please refer to the [trademark policy guidelines](http://www.openfoam.com/legal/trademark-policy.php).
+OpenCFD Ltd grants use of its OpenFOAM trademark by Third Parties on a licence basis. ESI Group and OpenFOAM Foundation Ltd are currently permitted to use the Name and agreed Domain Name. For information on trademark use, please refer to the [trademark policy guidelines](http://www.openfoam.com/legal/trademark-policy.php).
 
 Please [contact OpenCFD](http://www.openfoam.com/contact) if you have any questions on the use of the OpenFOAM trademark.
 

applications/solvers/basic/potentialFoam/overPotentialFoam/Make/files

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

applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -loverset

applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H

@@ -0,0 +1,9 @@
+const dictionary& potentialFlow
+(
+    mesh.solutionDict().subDict("potentialFlow")
+);
+
+const int nNonOrthCorr
+(
+    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
+);

applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H

@@ -0,0 +1,146 @@
+Info<< "Reading velocity field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Initialise the velocity internal field to zero
+U = dimensionedVector("0", U.dimensions(), Zero);
+
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    fvc::flux(U)
+);
+
+if (args.optionFound("initialiseUBCs"))
+{
+    U.correctBoundaryConditions();
+    phi = fvc::flux(U);
+}
+
+
+// Construct a pressure field
+// If it is available read it otherwise construct from the velocity BCs
+// converting fixed-value BCs to zero-gradient and vice versa.
+word pName("p");
+
+// Update name of the pressure field from the command-line option
+args.optionReadIfPresent("pName", pName);
+
+// Infer the pressure BCs from the velocity
+wordList pBCTypes
+(
+    U.boundaryField().size(),
+    fixedValueFvPatchScalarField::typeName
+);
+
+forAll(U.boundaryField(), patchi)
+{
+    if (U.boundaryField()[patchi].fixesValue())
+    {
+        pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
+    }
+}
+
+// Note that registerObject is false for the pressure field. The pressure
+// field in this solver doesn't have a physical value during the solution.
+// It shouldn't be looked up and used by sub models or boundary conditions.
+Info<< "Constructing pressure field " << pName << nl << endl;
+volScalarField p
+(
+    IOobject
+    (
+        pName,
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE,
+        false
+    ),
+    mesh,
+    dimensionedScalar(pName, sqr(dimVelocity), 0),
+    pBCTypes
+);
+
+// Infer the velocity potential BCs from the pressure
+wordList PhiBCTypes
+(
+    p.boundaryField().size(),
+    zeroGradientFvPatchScalarField::typeName
+);
+
+forAll(p.boundaryField(), patchi)
+{
+    if (p.boundaryField()[patchi].fixesValue())
+    {
+        PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
+    }
+}
+
+Info<< "Constructing velocity potential field Phi\n" << endl;
+volScalarField Phi
+(
+    IOobject
+    (
+        "Phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Phi", dimLength*dimVelocity, 0),
+    PhiBCTypes
+);
+
+label PhiRefCell = 0;
+scalar PhiRefValue = 0;
+setRefCell
+(
+    Phi,
+    potentialFlow.dict(),
+    PhiRefCell,
+    PhiRefValue
+);
+mesh.setFluxRequired(Phi.name());
+
+#include "createMRF.H"
+
+// Add overset specific interpolations
+{
+    dictionary oversetDict;
+    oversetDict.add("Phi", true);
+    oversetDict.add("U", true);
+
+    const_cast<dictionary&>
+    (
+        mesh.schemesDict()
+    ).add
+    (
+        "oversetInterpolationRequired",
+        oversetDict,
+        true
+    );
+}
+
+// Mask field for zeroing out contributions on hole cells
+#include "createCellMask.H"
+
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C

@@ -0,0 +1,260 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd
+     \\/     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
+    potentialFoam
+
+Group
+    grpBasicSolvers
+
+Description
+    Potential flow solver which solves for the velocity potential, to
+    calculate the flux-field, from which the velocity field is obtained by
+    reconstructing the flux.
+
+    \heading Solver details
+    The potential flow solution is typically employed to generate initial fields
+    for full Navier-Stokes codes.  The flow is evolved using the equation:
+
+    \f[
+        \laplacian \Phi = \div(\vec{U})
+    \f]
+
+    Where:
+    \vartable
+        \Phi      | Velocity potential [m2/s]
+        \vec{U}   | Velocity [m/s]
+    \endvartable
+
+    The corresponding pressure field could be calculated from the divergence
+    of the Euler equation:
+
+    \f[
+        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+
+    but this generates excessive pressure variation in regions of large
+    velocity gradient normal to the flow direction.  A better option is to
+    calculate the pressure field corresponding to velocity variation along the
+    stream-lines:
+
+    \f[
+        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+    where the flow direction tensor \f$\vec{F}\f$ is obtained from
+    \f[
+        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
+    \f]
+
+    \heading Required fields
+    \plaintable
+        U         | Velocity [m/s]
+    \endplaintable
+
+    \heading Optional fields
+    \plaintable
+        p         | Kinematic pressure [m2/s2]
+        Phi       | Velocity potential [m2/s]
+                  | Generated from p (if present) or U if not present
+    \endplaintable
+
+    \heading Options
+    \plaintable
+        -writep   | write the Euler pressure
+        -writePhi | Write the final velocity potential
+        -initialiseUBCs | Update the velocity boundaries before solving for Phi
+    \endplaintable
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pisoControl.H"
+#include "dynamicFvMesh.H"
+#include "cellCellStencilObject.H"
+#include "localMin.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "pName",
+        "pName",
+        "Name of the pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "initialiseUBCs",
+        "Initialise U boundary conditions"
+    );
+
+    argList::addBoolOption
+    (
+        "writePhi",
+        "Write the final velocity potential field"
+    );
+
+    argList::addBoolOption
+    (
+        "writep",
+        "Calculate and write the Euler pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "withFunctionObjects",
+        "execute functionObjects"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedDynamicFvMesh.H"
+
+    pisoControl potentialFlow(mesh, "potentialFlow");
+
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< nl << "Calculating potential flow" << endl;
+
+    mesh.update();
+
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+
+    // Since solver contains no time loop it would never execute
+    // function objects so do it ourselves
+    runTime.functionObjects().start();
+
+    MRF.makeRelative(phi);
+    adjustPhi(phi, U, p);
+
+    // Non-orthogonal velocity potential corrector loop
+    while (potentialFlow.correct())
+    {
+        phi = fvc::flux(U);
+
+        while (potentialFlow.correctNonOrthogonal())
+        {
+            fvScalarMatrix PhiEqn
+            (
+                fvm::laplacian(faceMask, Phi)
+            ==
+                fvc::div(phi)
+            );
+
+            PhiEqn.setReference(PhiRefCell, PhiRefValue);
+            PhiEqn.solve();
+
+            if (potentialFlow.finalNonOrthogonalIter())
+            {
+                phi -= PhiEqn.flux();
+            }
+        }
+
+        MRF.makeAbsolute(phi);
+
+        Info<< "Continuity error = "
+            << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
+            << endl;
+
+        U = fvc::reconstruct(phi);
+        U.correctBoundaryConditions();
+
+        Info<< "Interpolated velocity error = "
+            << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
+            << endl;
+    }
+
+    // Write U and phi
+    U.write();
+    phi.write();
+
+    // Optionally write Phi
+    if (args.optionFound("writePhi"))
+    {
+        Phi.write();
+    }
+
+    // Calculate the pressure field from the Euler equation
+    if (args.optionFound("writep"))
+    {
+        Info<< nl << "Calculating approximate pressure field" << endl;
+
+        label pRefCell = 0;
+        scalar pRefValue = 0.0;
+        setRefCell
+        (
+            p,
+            potentialFlow.dict(),
+            pRefCell,
+            pRefValue
+        );
+
+        // Calculate the flow-direction filter tensor
+        volScalarField magSqrU(magSqr(U));
+        volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
+
+        // Calculate the divergence of the flow-direction filtered div(U*U)
+        // Filtering with the flow-direction generates a more reasonable
+        // pressure distribution in regions of high velocity gradient in the
+        // direction of the flow
+        volScalarField divDivUU
+        (
+            fvc::div
+            (
+                F & fvc::div(phi, U),
+                "div(div(phi,U))"
+            )
+        );
+
+        // Solve a Poisson equation for the approximate pressure
+        while (potentialFlow.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::laplacian(p) + divDivUU
+            );
+
+            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.solve();
+        }
+
+        p.write();
+    }
+
+    runTime.functionObjects().end();
+
+    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C

@@ -40,15 +40,14 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
-        dictionaryConstructorTablePtr_->find(modelType);
 
     if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown PDRDragModel type "
             << modelType << nl << nl
-            << "Valid  PDRDragModels are : " << endl
+            << "Valid PDRDragModel types :" << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C

@@ -39,15 +39,14 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
-        dictionaryConstructorTablePtr_->find(modelType);
 
     if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiEqModel type "
             << modelType << nl << nl
-            << "Valid XiEqModels are : " << endl
+            << "Valid XiEqModel types :" << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C

@@ -39,15 +39,14 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
-        dictionaryConstructorTablePtr_->find(modelType);
 
     if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiGModel type "
             << modelType << nl << nl
-            << "Valid XiGModels are : " << endl
+            << "Valid XiGModel types :" << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C

@@ -42,15 +42,14 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
-        dictionaryConstructorTablePtr_->find(modelType);
 
     if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiModel type "
             << modelType << nl << nl
-            << "Valid XiModels are : " << endl
+            << "Valid XiModel types :" << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

applications/solvers/combustion/chemFoam/chemFoam.C

@@ -42,7 +42,7 @@ Description
 #include "OFstream.H"
 #include "thermoPhysicsTypes.H"
 #include "basicMultiComponentMixture.H"
-#include "cellModeller.H"
+#include "cellModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/combustion/chemFoam/createSingleCellMesh.H

@@ -13,8 +13,7 @@ points[5] = vector(1, 0, 1);
 points[6] = vector(1, 1, 1);
 points[7] = vector(0, 1, 1);
 
-const cellModel& hexa = *(cellModeller::lookup("hex"));
+faceList faces = cellModel::ref(cellModel::HEX).modelFaces();
-faceList faces = hexa.modelFaces();
 
 fvMesh mesh
 (
@@ -25,7 +24,7 @@ fvMesh mesh
         runTime,
         IOobject::READ_IF_PRESENT
     ),
-    xferMove<Field<vector>>(points),
+    xferMove<pointField>(points),
     faces.xfer(),
     owner.xfer(),
     neighbour.xfer()

applications/solvers/combustion/fireFoam/createFieldRefs.H

@@ -1,4 +1,4 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-filmModelType& surfaceFilm = tsurfaceFilm();
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
 const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H

@@ -1,5 +1,6 @@
 Info<< "\nConstructing surface film model" << endl;
 
-typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
-
+(
-autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
+    regionModels::surfaceFilmModel::New(mesh, g)
+);

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -209,16 +209,14 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
 
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    os.writeEntryIfDifferent<word>("U", "U", UName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
-    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
+    os.writeEntryIfDifferent<word>("mu", "thermo:mu", muName_);
 
-    os.writeKeyword("accommodationCoeff")
+    os.writeEntry("accommodationCoeff", accommodationCoeff_);
-        << accommodationCoeff_ << token::END_STATEMENT << nl;
     Twall_.writeEntry("Twall", os);
-    os.writeKeyword("gamma")
+    os.writeEntry("gamma", gamma_);
-        << gamma_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H

@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef smoluchowskiJumpTFvPatchScalarFields_H
+#ifndef smoluchowskiJumpTFvPatchScalarField_H
-#define smoluchowskiJumpTFvPatchScalarFields_H
+#define smoluchowskiJumpTFvPatchScalarField_H
 
 #include "mixedFvPatchFields.H"
 
@@ -43,7 +43,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                 Class smoluchowskiJumpTFvPatch Declaration
+            Class smoluchowskiJumpTFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class smoluchowskiJumpTFvPatchScalarField
@@ -74,6 +74,7 @@ class smoluchowskiJumpTFvPatchScalarField
         //- Heat capacity ratio (default 1.4)
         scalar gamma_;
 
+
 public:
 
     //- Runtime type information

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C

@@ -200,18 +200,16 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()
 void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
-    writeEntryIfDifferent<word>(os, "T", "T", TName_);
+    os.writeEntryIfDifferent<word>("T", "T", TName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
-    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
+    os.writeEntryIfDifferent<word>("mu", "thermo:mu", muName_);
-    writeEntryIfDifferent<word>(os, "tauMC", "tauMC", tauMCName_);
+    os.writeEntryIfDifferent<word>("tauMC", "tauMC", tauMCName_);
 
-    os.writeKeyword("accommodationCoeff")
+    os.writeEntry("accommodationCoeff", accommodationCoeff_);
-        << accommodationCoeff_ << token::END_STATEMENT << nl;
     Uwall_.writeEntry("Uwall", os);
-    os.writeKeyword("thermalCreep")
+    os.writeEntry("thermalCreep", thermalCreep_);
-        << thermalCreep_ << token::END_STATEMENT << nl;
+    os.writeEntry("curvature", curvature_);
-    os.writeKeyword("curvature") << curvature_ << token::END_STATEMENT << nl;
 
     refValue().writeEntry("refValue", os);
     valueFraction().writeEntry("valueFraction", os);

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C

@@ -117,8 +117,8 @@ void Foam::fixedRhoFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
 
-    writeEntryIfDifferent<word>(os, "p", "p", this->pName_);
+    os.writeEntryIfDifferent<word>("p", "p", pName_);
-    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
     writeEntry("value", os);
 }
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -113,6 +113,9 @@ int main(int argc, char *argv[])
             phiv_pos -= mesh.phi();
             phiv_neg -= mesh.phi();
         }
+        // Note: extracted out the orientation so becomes unoriented
+        phiv_pos.setOriented(false);
+        phiv_neg.setOriented(false);
 
         volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
@@ -120,14 +123,11 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (
@@ -168,11 +168,12 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
-        surfaceVectorField phiUp
+        surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
-        (
+        // Note: reassembled orientation from the pos and neg parts so becomes
-            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
+        // oriented
-          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        phiU.setOriented(true);
-        );
+
+        surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
 
         surfaceScalarField phiEp
         (
@@ -185,7 +186,10 @@ int main(int argc, char *argv[])
         // Make flux for pressure-work absolute
         if (mesh.moving())
         {
-            phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
+            surfaceScalarField phia(a_pos*p_pos + a_neg*p_neg);
+            phia.setOriented(true);
+
+            phiEp += mesh.phi()*phia;
         }
 
         volScalarField muEff("muEff", turbulence->muEff());

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -93,7 +93,10 @@ int main(int argc, char *argv[])
         surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
 
         surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
+        // Note: extracted out the orientation so becomes unoriented
+        phiv_pos.setOriented(false);
         surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
+        phiv_neg.setOriented(false);
 
         volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
@@ -101,20 +104,19 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
 
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (
             "ap",
             max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
         );
+
         surfaceScalarField am
         (
             "am",
@@ -163,11 +165,12 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
-        surfaceVectorField phiUp
+        surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
-        (
+        // Note: reassembled orientation from the pos and neg parts so becomes
-            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
+        // oriented
-          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        phiU.setOriented(true);
-        );
+
+        surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
 
         surfaceScalarField phiEp
         (

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H

@@ -61,18 +61,6 @@ dimensionedScalar rhoMin
     )
 );
 
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
 mesh.setFluxRequired(p.name());
 
 Info<< "Creating field dpdt\n" << endl;
@@ -115,3 +103,15 @@ volScalarField K("K", 0.5*magSqr(U));
 
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H

@@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
+volVectorField HbyA("HbyA", U);
-//mesh.interpolate(HbyA);
+HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude
+
+EXE_LIBS = \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfvOptions \
+    -loverset \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/UEqn.H

@@ -0,0 +1,23 @@
+    // Solve the Momentum equation
+    MRF.correctBoundaryVelocity(U);
+
+    tmp<fvVectorMatrix> tUEqn
+    (
+        fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        fvOptions(rho, U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (simple.momentumPredictor())
+    {
+        solve(UEqn == -cellMask*fvc::grad(p));
+    }
+
+    fvOptions.correct(U);

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFieldRefs.H

@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H

@@ -0,0 +1,91 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+autoPtr<fluidThermo> pThermo
+(
+    fluidThermo::New(mesh)
+);
+fluidThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+pressureControl pressureControl(p, rho, simple.dict());
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+#include "createMRF.H"
+
+//- Overset specific
+
+// Add solver-specific interpolations
+{
+    dictionary oversetDict;
+    oversetDict.add("U", true);
+    oversetDict.add("p", true);
+    oversetDict.add("HbyA", true);
+    oversetDict.add("grad(p)", true);
+    oversetDict.add("rho", true);
+
+    const_cast<dictionary&>
+    (
+        mesh.schemesDict()
+    ).add
+    (
+        "oversetInterpolationRequired",
+        oversetDict,
+        true
+    );
+}
+
+// Mask field for zeroing out contributions on hole cells
+#include "createCellMask.H"
+
+#include "createInterpolatedCells.H"
+
+bool adjustFringe
+(
+    simple.dict().lookupOrDefault("oversetAdjustPhi", false)
+);

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H

@@ -0,0 +1,22 @@
+    Info<< "Create dynamic mesh for time = "
+        << runTime.timeName() << nl << endl;
+
+    autoPtr<dynamicFvMesh> meshPtr
+    (
+        dynamicFvMesh::New
+        (
+            IOobject
+            (
+                dynamicFvMesh::defaultRegion,
+                runTime.timeName(),
+                runTime,
+                IOobject::MUST_READ
+            )
+        )
+    );
+
+    dynamicFvMesh& mesh = meshPtr();
+
+    // Calculate initial mesh-to-mesh mapping. Note that this should be
+    // done under the hood, e.g. as a MeshObject
+    mesh.update();

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/pEqn.H

@@ -0,0 +1,123 @@
+{
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
+    volVectorField HbyA("HbyA", U);
+    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
+    tUEqn.clear();
+
+    bool closedVolume = false;
+
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+    if (simple.transonic())
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+        );
+
+        phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+        while (simple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::div(phiHbyA)
+              + fvm::div(phid, p)
+              - fvm::laplacian(rhorAUf, p)
+            ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            // Relax the pressure equation to ensure diagonal-dominance
+            pEqn.relax();
+
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        closedVolume = adjustPhi(phiHbyA, U, p);
+
+        if (adjustFringe)
+        {
+            oversetAdjustPhi(phiHbyA, U);
+        }
+
+        while (simple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::div(phiHbyA)
+              - fvm::laplacian(rhorAUf, p)
+            ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
+        }
+    }
+
+    #include "incompressible/continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    volVectorField gradP(fvc::grad(p));
+
+    U = HbyA - rAU*cellMask*gradP;
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+
+    bool pLimited = pressureControl.limit(p);
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+    }
+
+    if (pLimited || closedVolume)
+    {
+        p.correctBoundaryConditions();
+    }
+
+    rho = thermo.rho();
+
+    if (!simple.transonic())
+    {
+        rho.relax();
+    }
+}

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C

@@ -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) 2017 OpenCFD Ltd
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,43 +22,41 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    reactingParcelFilmFoam
+    overRhoSimpleFoam
 
 Group
-    grpLagrangianSolvers
+    grpCompressibleSolvers
 
 Description
-    Transient solver for compressible, turbulent flow with a reacting,
+    Overset steady-state solver for turbulent flow of compressible fluids.
-    multiphase particle cloud, and surface film modelling.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingCloud.H"
+#include "simpleControl.H"
-#include "surfaceFilmModel.H"
+#include "pressureControl.H"
-#include "psiCombustionModel.H"
-#include "radiationModel.H"
-#include "SLGThermo.H"
 #include "fvOptions.H"
-#include "pimpleControl.H"
+#include "cellCellStencilObject.H"
+#include "localMin.H"
+#include "oversetAdjustPhi.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    #define CREATE_MESH createUpdatedDynamicFvMesh.H
     #include "postProcess.H"
 
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createMesh.H"
+    #include "createUpdatedDynamicFvMesh.H"
     #include "createControl.H"
-    #include "createTimeControls.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createFvOptions.H"
-    #include "compressibleCourantNo.H"
-    #include "setInitialDeltaT.H"
     #include "initContinuityErrs.H"
 
     turbulence->validate();
@@ -67,46 +65,16 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (simple.loop())
     {
-        #include "readTimeControls.H"
-        #include "compressibleCourantNo.H"
-        #include "setMultiRegionDeltaT.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        parcels.evolve();
+        // Pressure-velocity SIMPLE corrector
-
-        surfaceFilm.evolve();
-
-        if (solvePrimaryRegion)
-        {
-            #include "rhoEqn.H"
-
-            // --- PIMPLE loop
-            while (pimple.loop())
-            {
         #include "UEqn.H"
-                #include "YEqn.H"
         #include "EEqn.H"
-
-                // --- Pressure corrector loop
-                while (pimple.correct())
-                {
         #include "pEqn.H"
-                }
 
-                if (pimple.turbCorr())
-                {
         turbulence->correct();
-                }
-            }
-
-            rho = thermo.rho();
-        }
 
         runTime.write();
 
@@ -115,7 +83,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/solvers/doc/solver.dox

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ The available solvers are grouped into the following categories:
   - \ref grpLagrangianSolvers
   - \ref grpMultiphaseSolvers
   - \ref grpStressAnalysisSolvers
+  - \ref grpFiniteAreaSolvers
 
 
 \*---------------------------------------------------------------------------*/

applications/solvers/doc/solversDoc.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,4 +34,10 @@ License
     This group contains moving mesh solvers solvers
 @}
 
+\defgroup grpFiniteAreaSolvers Finite area solvers
+@{
+    \ingroup grpSolvers
+    This group contains finite area solvers
+@}
+
 \*---------------------------------------------------------------------------*/

applications/solvers/finiteArea/liquidFilmFoam/Make/files

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

applications/solvers/finiteArea/liquidFilmFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/liquidFilmFoam/calcFrictionFactor.H

@@ -0,0 +1,7 @@
+{
+    // Stabilisation of friction factor calculation
+    // Friction factor is defined with standard gravity
+    frictionFactor.primitiveFieldRef() =
+        mag(2*9.81*sqr(manningField.primitiveField())/
+        pow(mag(h.primitiveField()) + 1e-7, 1.0/3.0));
+}

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

@@ -0,0 +1,13 @@
+{
+    scalar CoNumSigma = max
+    (
+        sqrt
+        (
+            2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
+            *aMesh.edgeInterpolation::deltaCoeffs()
+            /rhol
+        )
+    ).value()*runTime.deltaT().value();
+
+    Info<< "Max Capillary Courant Number = " << CoNumSigma << '\n' << endl;
+}

applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H

@@ -0,0 +1,158 @@
+    Info<< "Reading field h" << endl;
+    areaScalarField h
+    (
+        IOobject
+        (
+            "h",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+
+    Info<< "Reading field Us" << endl;
+    areaVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+
+    edgeScalarField phis
+    (
+        IOobject
+        (
+            "phis",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(Us) & aMesh.Le()
+    );
+
+
+    edgeScalarField phi2s
+    (
+        IOobject
+        (
+            "phi2s",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(h*Us) & aMesh.Le()
+    );
+
+
+    const areaVectorField& Ns = aMesh.faceAreaNormals();
+    areaVectorField Gs(g - Ns*(Ns & g));
+    areaScalarField Gn(mag(g - Gs));
+
+    // Mass source
+    areaScalarField Sm
+    (
+        IOobject
+        (
+            "Sm",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("Sm", dimLength/dimTime, 0)
+    );
+
+    // Mass sink
+    areaScalarField Sd
+    (
+        IOobject
+        (
+            "Sd",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("Sd", dimLength/dimTime, 0)
+    );
+
+    areaVectorField Ug
+    (
+        IOobject
+        (
+            "Sg",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedVector("Ug", dimVelocity, vector::zero)
+    );
+
+
+    // Surface pressure
+    areaScalarField ps
+    (
+        IOobject
+        (
+            "ps",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        rhol*Gn*h - sigma*fac::laplacian(h)
+    );
+
+    // Friction factor
+    areaScalarField dotProduct
+    (
+        aMesh.faceAreaNormals() & (g/mag(g))
+    );
+
+    Info<< "View factor: min = " << min(dotProduct.internalField())
+        << " max = " << max(dotProduct.internalField()) << endl;
+
+    areaScalarField manningField
+    (
+        IOobject
+        (
+            "manningField",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        aMesh
+    );
+
+    areaScalarField frictionFactor
+    (
+        IOobject
+        (
+            "frictionFactor",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar("one", dimless, 0.01)
+    );
+
+    aMesh.setFluxRequired("h");

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

@@ -0,0 +1,31 @@
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", dimVelocity, vector::zero)
+);
+
+
+volScalarField H
+(
+    IOobject
+    (
+        "H",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimLength, 0)
+);
+
+// Create volume-to surface mapping object
+volSurfaceMapping vsm(aMesh);

applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+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
+    liquidFilmFoam
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids in
+    liquid film formulation.
+
+Author
+    Zeljko Tukovic, FMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "loopControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFaMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readTransportProperties.H"
+    #include "createFaFields.H"
+    #include "createFvFields.H"
+    #include "createTimeControls.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readSolutionControls.H"
+        #include "readTimeControls.H"
+        #include "surfaceCourantNo.H"
+        #include "capillaryCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        while (iters.loop())
+        {
+            phi2s = fac::interpolate(h)*phis;
+
+            #include "calcFrictionFactor.H"
+
+            faVectorMatrix UsEqn
+            (
+                fam::ddt(h, Us)
+              + fam::div(phi2s, Us)
+              + fam::Sp(0.0125*frictionFactor*mag(Us), Us)
+             ==
+                Gs*h
+              - fam::Sp(Sd, Us)
+            );
+
+            UsEqn.relax();
+            solve(UsEqn == - fac::grad(ps*h)/rhol + ps*fac::grad(h)/rhol);
+
+            areaScalarField UsA(UsEqn.A());
+
+            Us = UsEqn.H()/UsA;
+            Us.correctBoundaryConditions();
+
+            phis =
+                (fac::interpolate(Us) & aMesh.Le())
+              - fac::interpolate(1.0/(rhol*UsA))*fac::lnGrad(ps*h)*aMesh.magLe()
+              + fac::interpolate(ps/(rhol*UsA))*fac::lnGrad(h)*aMesh.magLe();
+
+            faScalarMatrix hEqn
+            (
+                fam::ddt(h)
+              + fam::div(phis, h)
+             ==
+                Sm
+              - fam::Sp
+                (
+                    Sd/(h + dimensionedScalar("small", dimLength, SMALL)),
+                    h
+                )
+            );
+
+            hEqn.relax();
+            hEqn.solve();
+
+            phi2s = hEqn.flux();
+
+            // Bound h
+            h.primitiveFieldRef() = max
+            (
+                max
+                (
+                    h.primitiveField(),
+                    fac::average(max(h, h0))().primitiveField()
+                   *pos(h0.value() - h.primitiveField())
+                ),
+                h0.value()
+            );
+
+            ps = rhol*Gn*h - sigma*fac::laplacian(h);
+            ps.correctBoundaryConditions();
+
+            Us -= (1.0/(rhol*UsA))*fac::grad(ps*h)
+                - (ps/(rhol*UsA))*fac::grad(h);
+            Us.correctBoundaryConditions();
+        }
+
+        if (runTime.outputTime())
+        {
+            vsm.mapToVolume(h, H.boundaryFieldRef());
+            vsm.mapToVolume(Us, U.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H

@@ -0,0 +1 @@
+loopControl iters(runTime, aMesh.solutionDict(), "solution");

applications/solvers/finiteArea/liquidFilmFoam/readTransportProperties.H

@@ -0,0 +1,41 @@
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+dimensionedScalar mug
+(
+    transportProperties.lookup("mug")
+);
+
+dimensionedScalar mul
+(
+    transportProperties.lookup("mul")
+);
+
+dimensionedScalar sigma
+(
+    transportProperties.lookup("sigma")
+);
+
+dimensionedScalar rhol
+(
+    transportProperties.lookup("rhol")
+);
+
+dimensionedScalar rhog
+(
+    transportProperties.lookup("rhog")
+);
+
+dimensionedScalar h0
+(
+    transportProperties.lookup("h0")
+);

applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+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/>.
+
+Global
+    surfaceCourantNo
+
+Author
+    Hrvoje Jasak, Wikki Ltd.
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers for the
+    Finite Area method.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum = 0.0;
+scalar meanCoNum = 0.0;
+scalar velMag = 0.0;
+
+if (aMesh.nInternalEdges())
+{
+    edgeScalarField SfUfbyDelta
+    (
+        aMesh.edgeInterpolation::deltaCoeffs()*mag(phis)
+    );
+
+    CoNum = max(SfUfbyDelta/aMesh.magLe())
+        .value()*runTime.deltaT().value();
+
+    meanCoNum = (sum(SfUfbyDelta)/sum(aMesh.magLe()))
+        .value()*runTime.deltaT().value();
+
+    velMag = max(mag(phis)/aMesh.magLe()).value();
+}
+
+Info<< "Courant Number mean: " << meanCoNum
+    << " max: " << CoNum
+    << " velocity magnitude: " << velMag << endl;
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/sphereSurfactantFoam/Make/files

@@ -0,0 +1,3 @@
+surfactantFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/sphereSurfactantFoam

applications/solvers/finiteArea/sphereSurfactantFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H

@@ -0,0 +1,78 @@
+Info<< "Reading field Cs" << endl;
+areaScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    aMesh
+);
+
+dimensioned<scalar> Cs0
+(
+    "Cs0",
+    dimensionSet(1, -2, 0, 0, 0, 0, 0),
+    1.0
+);
+
+const areaVectorField& R = aMesh.areaCentres();
+
+Cs = Cs0*(1.0 + R.component(vector::X)/mag(R));
+
+
+dimensioned<scalar> Ds
+(
+    "Ds",
+    dimensionSet(0, 2, -1, 0, 0, 0, 0),
+    1.0
+);
+
+
+areaVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    aMesh,
+    dimensioned<vector>("Us", dimVelocity, vector::zero)
+);
+
+dimensioned<scalar> Uinf("Uinf", dimVelocity, 1.0);
+
+forAll (Us, faceI)
+{
+    Us[faceI].x() =
+        Uinf.value()*(0.25*(3.0 + sqr(R[faceI].x()/mag(R[faceI]))) - 1.0);
+
+    Us[faceI].y() =
+        Uinf.value()*0.25*R[faceI].x()*R[faceI].y()/sqr(mag(R[faceI]));
+
+    Us[faceI].z() =
+        Uinf.value()*0.25*R[faceI].x()*R[faceI].z()/sqr(mag(R[faceI]));
+}
+
+Us -= aMesh.faceAreaNormals()*(aMesh.faceAreaNormals() & Us);
+
+
+edgeScalarField phis
+(
+    IOobject
+    (
+        "phis",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearEdgeInterpolate(Us) & aMesh.Le()
+);
+

applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

@@ -0,0 +1,36 @@
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
+
+    volScalarField Cvf
+    (
+        IOobject
+        (
+            "Cvf",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless/dimLength, 0)
+    );
+
+    vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+    Cvf.write();
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimVelocity, vector::zero)
+    );
+
+    vsm.mapToVolume(Us, U.boundaryFieldRef());
+    U.write();

applications/solvers/finiteArea/sphereSurfactantFoam/surfactantFoam.C

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+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
+    surfactantFoam for sphere transport
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Passive scalar transport equation solver on a sphere
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    #include "createFaMesh.H"
+    #include "createFaFields.H"
+    #include "createVolFields.H"
+
+    Info<< "Total mass of surfactant: "
+        << sum(Cs.internalField()*aMesh.S()) << endl;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.value() << endl;
+
+        faScalarMatrix CsEqn
+        (
+            fam::ddt(Cs)
+          + fam::div(phis, Cs)
+          - fam::laplacian(Ds, Cs)
+        );
+
+        CsEqn.solve();
+
+        if (runTime.writeTime())
+        {
+            vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "Total mass of surfactant: "
+            << sum(Cs.internalField()*aMesh.S()) << endl;
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/finiteArea/surfactantFoam/Make/files

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

applications/solvers/finiteArea/surfactantFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/cfdTools/general/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/finiteArea/surfactantFoam/createFaFields.H

@@ -0,0 +1,63 @@
+Info<< "Reading field Cs" << endl;
+areaScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    aMesh
+);
+
+Info<< "Reading transportProperties\n" << endl;
+
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+
+Info<< "Reading diffusivity D\n" << endl;
+
+dimensionedScalar Ds
+(
+    transportProperties.lookup("Ds")
+);
+
+areaVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    aMesh
+);
+
+
+edgeScalarField phis
+(
+    IOobject
+    (
+        "phis",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearEdgeInterpolate(Us) & aMesh.Le()
+);
+

applications/solvers/finiteArea/surfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

@@ -0,0 +1,36 @@
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
+
+    volScalarField Cvf
+    (
+        IOobject
+        (
+            "Cvf",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless/dimLength, 0)
+    );
+
+    vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+    Cvf.write();
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimVelocity, vector::zero)
+    );
+
+    vsm.mapToVolume(Us, U.boundaryFieldRef());
+    U.write();

applications/solvers/finiteArea/surfactantFoam/surfactantFoam.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2016-2017 Wikki Ltd
+-------------------------------------------------------------------------------
+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
+    surfactantFoam
+
+Group
+    grpFiniteAreaSolvers
+
+Description
+    Passive scalar transport equation solver.
+
+    \heading Solver details
+    The equation is given by:
+
+    \f[
+        \ddt{Cs} + \div \left(\vec{U} Cs\right) - \div \left(D_T \grad Cs \right)
+        = 0
+    \f]
+
+    Where:
+    \vartable
+        Cs      | Passive scalar
+        Ds      | Diffusion coefficient
+    \endvartable
+
+    \heading Required fields
+    \plaintable
+        Cs      | Passive scalar
+        Us      | Velocity [m/s]
+    \endplaintable
+
+Author
+    Zeljko Tukovic, FMENA
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFaMesh.H"
+    #include "createFaFields.H"
+    #include "createVolFields.H"
+
+    Info<< "Total mass of surfactant: "
+        << sum(Cs.internalField()*aMesh.S()) << endl;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.value() << endl;
+
+        faScalarMatrix CsEqn
+        (
+            fam::ddt(Cs)
+          + fam::div(phis, Cs)
+          - fam::laplacian(Ds, Cs)
+        );
+
+        CsEqn.solve();
+
+        if (runTime.writeTime())
+        {
+            vsm.mapToVolume(Cs, Cvf.boundaryFieldRef());
+
+            runTime.write();
+        }
+
+        Info<< "Total mass of surfactant: "
+            << sum(Cs.internalField()*aMesh.S()) << endl;
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,7 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
+#include "loopControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,11 +90,10 @@ int main(int argc, char *argv[])
             }
         }
 
-
         // --- PIMPLE loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (int oCorr=0; oCorr<nOuterCorr; ++oCorr)
         {
-            bool finalIter = oCorr == nOuterCorr-1;
+            const bool finalIter = (oCorr == nOuterCorr-1);
 
             forAll(fluidRegions, i)
             {
@@ -113,6 +113,35 @@ int main(int argc, char *argv[])
                 #include "solveSolid.H"
             }
 
+            // Additional loops for energy solution only
+            if (!oCorr && nOuterCorr > 1)
+            {
+                loopControl looping(runTime, pimple, "energyCoupling");
+
+                while (looping.loop())
+                {
+                    Info<< nl << looping << nl;
+
+                    forAll(fluidRegions, i)
+                    {
+                        Info<< "\nSolving for fluid region "
+                            << fluidRegions[i].name() << endl;
+                       #include "setRegionFluidFields.H"
+                       #include "readFluidMultiRegionPIMPLEControls.H"
+                       frozenFlow = true;
+                       #include "solveFluid.H"
+                    }
+
+                    forAll(solidRegions, i)
+                    {
+                        Info<< "\nSolving for solid region "
+                            << solidRegions[i].name() << endl;
+                        #include "setRegionSolidFields.H"
+                        #include "readSolidMultiRegionPIMPLEControls.H"
+                        #include "solveSolid.H"
+                    }
+                }
+            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
+#include "loopControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +58,6 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
-
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -80,6 +80,35 @@ int main(int argc, char *argv[])
             #include "solveSolid.H"
         }
 
+        // Additional loops for energy solution only
+        {
+            loopControl looping(runTime, "SIMPLE", "energyCoupling");
+
+            while (looping.loop())
+            {
+                Info<< nl << looping << nl;
+
+                forAll(fluidRegions, i)
+                {
+                    Info<< "\nSolving for fluid region "
+                        << fluidRegions[i].name() << endl;
+                   #include "setRegionFluidFields.H"
+                   #include "readFluidMultiRegionSIMPLEControls.H"
+                   frozenFlow = true;
+                   #include "solveFluid.H"
+                }
+
+                forAll(solidRegions, i)
+                {
+                    Info<< "\nSolving for solid region "
+                        << solidRegions[i].name() << endl;
+                    #include "setRegionSolidFields.H"
+                    #include "readSolidMultiRegionSIMPLEControls.H"
+                    #include "solveSolid.H"
+                }
+            }
+        }
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -12,6 +12,14 @@
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
 
+    volScalarField& p_rgh = p_rghFluid[i];
+
+    const dimensionedVector& g = gFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];
+
+    radiation::radiationModel& rad = radiation[i];
+
     IOMRFZoneList& MRF = MRFfluid[i];
     fv::options& fvOptions = fluidFvOptions[i];
 
@@ -22,14 +30,7 @@
         initialMassFluid[i]
     );
 
-    radiation::radiationModel& rad = radiation[i];
+    bool frozenFlow = frozenFlowFluid[i];
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
-    const bool frozenFlow = frozenFlowFluid[i];
-
-    volScalarField& p_rgh = p_rghFluid[i];
-
-    const dimensionedVector& g = gFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H

@@ -1,5 +1,5 @@
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
     {
         fvScalarMatrix hEqn
         (
@@ -20,9 +20,9 @@
 
         fvOptions.correct(h);
     }
-}
 
     thermo.correct();
 
     Info<< "Min/max T:" << min(thermo.T()).value() << ' '
         << max(thermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -19,8 +19,8 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
-List<label> refCellFluid(fluidRegions.size());
+List<label> pRefCellFluid(fluidRegions.size());
-List<scalar> refValueFluid(fluidRegions.size());
+List<scalar> pRefValueFluid(fluidRegions.size());
 
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
@@ -252,8 +252,8 @@ forAll(fluidRegions, i)
 
     turbulence[i].validate();
 
-    refCellFluid[i] = 0;
+    pRefCellFluid[i] = 0;
-    refValueFluid[i] = 0.0;
+    pRefValueFluid[i] = 0.0;
 
     if (p_rghFluid[i].needReference())
     {
@@ -262,8 +262,8 @@ forAll(fluidRegions, i)
             thermoFluid[i].p(),
             p_rghFluid[i],
             pimpleDict,
-            refCellFluid[i],
+            pRefCellFluid[i],
-            refValueFluid[i]
+            pRefValueFluid[i]
         );
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -32,7 +32,8 @@
         initialMassFluid[i]
     );
 
-    const bool frozenFlow = frozenFlowFluid[i];
+    bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = pRefCellFluid[i];
+    const scalar pRefValue = pRefValueFluid[i];
 
-    const label pRefCell = refCellFluid[i];
-    const scalar pRefValue = refValueFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H

@@ -1,10 +1,10 @@
-    const wordList solidsNames(rp["solid"]);
+    const wordList solidNames(rp["solid"]);
 
-    PtrList<fvMesh> solidRegions(solidsNames.size());
+    PtrList<fvMesh> solidRegions(solidNames.size());
 
-    forAll(solidsNames, i)
+    forAll(solidNames, i)
     {
-        Info<< "Create solid mesh for region " << solidsNames[i]
+        Info<< "Create solid mesh for region " << solidNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         solidRegions.set
@@ -14,7 +14,7 @@
             (
                 IOobject
                 (
-                    solidsNames[i],
+                    solidNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H

@@ -27,8 +27,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef solidRegionDiff_H
+#ifndef solidRegionDiffNo_H
-#define solidRegionDiff_H
+#define solidRegionDiffNo_H
 
 #include "fvMesh.H"
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -4,7 +4,7 @@ if (finalIter)
 }
 
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth)
     {
         fvScalarMatrix hEqn
         (
@@ -26,12 +26,12 @@ if (finalIter)
 
         fvOptions.correct(h);
     }
-}
 
     thermo.correct();
 
     Info<< "Min/max T:" << min(thermo.T()).value() << ' '
         << max(thermo.T()).value() << endl;
+}
 
 if (finalIter)
 {

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H

@@ -20,7 +20,7 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
+HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
 //mesh.interpolate(HbyA);
 if (massFluxInterpolation)

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pimpleDyMFoam.C

@@ -54,10 +54,8 @@ Description
 
 int main(int argc, char *argv[])
 {
-    argList::addNote
+    #include "postProcess.H"
-    (
+
-        "Experimental version of pimpleDyMFoam with support for overset meshes"
-    );
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"

applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H

@@ -16,6 +16,9 @@
 
     fvOptions.constrain(UEqn);
 
+    if (simple.momentumPredictor())
+    {
         solve(UEqn == -cellMask*fvc::grad(p));
+    }
 
     fvOptions.correct(U);

applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H

@@ -5,7 +5,7 @@
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
+    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
     //mesh.interpolate(HbyA);
     if (massFluxInterpolation)

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

@@ -77,6 +77,9 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        // Store the particle positions
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         // Calculate absolute flux from the mapped surface velocity

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C

@@ -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
@@ -68,6 +68,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         U.correctBoundaryConditions();

applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files

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

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -1,34 +0,0 @@
-    MRF.correctBoundaryVelocity(U);
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + MRF.DDt(rho, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        parcels.SU(U)
-      + fvOptions(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-          ==
-            fvc::reconstruct
-            (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                )*mesh.magSf()
-            )
-        );
-
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }

applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H

@@ -1,9 +0,0 @@
-Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
-(
-    "reactingCloud1",
-    rho,
-    U,
-    g,
-    slgThermo
-);

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -1,141 +0,0 @@
-Info<< "Creating combustion model\n" << endl;
-
-autoPtr<combustionModels::psiCombustionModel> combustion
-(
-    combustionModels::psiCombustionModel::New(mesh)
-);
-
-psiReactionThermo& thermo = combustion->thermo();
-thermo.validate(args.executable(), "h", "e");
-
-SLGThermo slgThermo(mesh, thermo);
-
-basicSpecieMixture& composition = thermo.composition();
-PtrList<volScalarField>& Y = composition.Y();
-
-const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
-
-Info<< "Creating field rho\n" << endl;
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "\nReading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-// Set the turbulence into the combustion model
-combustion->setTurbulence(turbulence());
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh;
-
-mesh.setFluxRequired(p_rgh.name());
-
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-forAll(Y, i)
-{
-    fields.add(Y[i]);
-}
-fields.add(thermo.he());
-
-IOdictionary additionalControlsDict
-(
-    IOobject
-    (
-        "additionalControls",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
-
-Switch solvePrimaryRegion
-(
-    additionalControlsDict.lookup("solvePrimaryRegion")
-);
-
-volScalarField Qdot
-(
-    IOobject
-    (
-        "Qdot",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createMRF.H"
-#include "createClouds.H"
-#include "createRadiationModel.H"
-#include "createSurfaceFilmModel.H"

applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H

@@ -1,5 +0,0 @@
-Info<< "\nConstructing surface film model" << endl;
-
-typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
-
-autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,59 +0,0 @@
-rho = thermo.rho();
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
-  + phig
-);
-
-MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-while (pimple.correctNonOrthogonal())
-{
-    fvScalarMatrix p_rghEqn
-    (
-        fvc::ddt(psi, rho)*gh
-      + fvc::div(phiHbyA)
-      + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(rhorAUf, p_rgh)
-     ==
-        parcels.Srho()
-      + surfaceFilm.Srho()
-      + fvOptions(psi, p_rgh, rho.name())
-    );
-
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-    }
-}
-
-p = p_rgh + rho*gh;
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-K = 0.5*magSqr(U);
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}

applications/solvers/lagrangian/reactingParcelFoam/EEqn.H

@@ -19,6 +19,7 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
+      + surfaceFilm.Sh()
       + radiation->Sh(thermo, he)
       + Qdot
       + fvOptions(rho, he)
@@ -35,6 +36,6 @@
     thermo.correct();
     radiation->correct();
 
-    Info<< "T gas min/max   " << min(T).value() << ", "
+    Info<< "T gas min/max   = " << min(T).value() << ", "
         << max(T).value() << endl;
 }

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -1,12 +1,11 @@
 EXE_INC = \
     -I. \
+    -I../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -17,33 +16,33 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
-
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
+    -lsampling \
     -lmeshTools \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -llagrangianTurbulence \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \
-    -lradiationModels \
-    -lODE \
     -lregionModels \
+    -lradiationModels \
     -lsurfaceFilmModels \
-    -lcombustionModels \
+    -lsurfaceFilmDerivedFvPatchFields \
-    -lfvOptions \
+    -llagrangian \
-    -lsampling
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lODE \
+    -lcombustionModels

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -6,8 +6,7 @@
       + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
-        rho()*g
+        parcels.SU(U)
-      + parcels.SU(U)
       + fvOptions(rho, U)
     );
 
@@ -17,7 +16,18 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p));
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
 
         fvOptions.correct(U);
         K = 0.5*magSqr(U);

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -9,6 +9,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 
+
 {
     combustion->correct();
     Qdot = combustion->Qdot();
@@ -24,11 +25,12 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
               ==
                 parcels.SYi(i, Yi)
-              + combustion->R(Yi)
               + fvOptions(rho, Yi)
+              + combustion->R(Yi)
+              + surfaceFilm.Srho(i)
             );
 
             YEqn.relax();

applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H

@@ -1,3 +1,5 @@
+const label inertIndex(composition.species()[inertSpecie]);
+
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-const label inertIndex(composition.species()[inertSpecie]);
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,7 +1,5 @@
 #include "createRDeltaT.H"
 
-#include "readGravitationalAcceleration.H"
-
 Info<< "Creating combustion model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> combustion
@@ -26,8 +24,7 @@ if (!composition.species().found(inertSpecie))
         << exit(FatalIOError);
 }
 
-volScalarField& p = thermo.p();
+Info<< "Creating field rho\n" << endl;
-
 volScalarField rho
 (
     IOobject
@@ -41,6 +38,8 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField& p = thermo.p();
+
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
@@ -57,30 +56,6 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -96,6 +71,31 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
+mesh.setFluxRequired(p_rgh.name());
+
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -105,6 +105,11 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
+Switch solvePrimaryRegion
+(
+    pimple.dict().lookupOrDefault<Switch>("solvePrimaryRegion", true)
+);
+
 volScalarField Qdot
 (
     IOobject
@@ -126,3 +131,4 @@ volScalarField Qdot
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"
+#include "createSurfaceFilmModel.H"

applications/solvers/lagrangian/reactingParcelFoam/createSurfaceFilmModel.H

@@ -0,0 +1,6 @@
+Info<< "\nConstructing surface film model" << endl;
+
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
+(
+    regionModels::surfaceFilmModel::New(mesh, g)
+);

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -1,4 +1,7 @@
+if (!pimple.SIMPLErho())
+{
     rho = thermo.rho();
+}
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
 // pressure solution
@@ -7,6 +10,9 @@ const volScalarField psip0(psi*p);
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -14,58 +20,68 @@ surfaceScalarField phiHbyA
         fvc::flux(rho*HbyA)
       + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
+  + phig
 );
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the pressure BCs to ensure flux consistency
-constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-fvScalarMatrix pDDtEqn
+fvScalarMatrix p_rghDDtEqn
 (
-    fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
   + fvc::div(phiHbyA)
  ==
     parcels.Srho()
-  + fvOptions(psi, p, rho.name())
+  + surfaceFilm.Srho()
+  + fvOptions(psi, p_rgh, rho.name())
 );
 
 while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix pEqn
+    fvScalarMatrix p_rghEqn
     (
-        pDDtEqn
+        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p)
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi = phiHbyA + pEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
 }
 
-p.relax();
+p = p_rgh + rho*gh;
 
 // Thermodynamic density update
 thermo.correctRho(psi*p - psip0);
 
-#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+#include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U = HbyA - rAU*fvc::grad(p);
+if (pressureControl.limit(p))
-U.correctBoundaryConditions();
+{
-fvOptions.correct(U);
+    p.correctBoundaryConditions();
-K = 0.5*magSqr(U);
+    rho = thermo.rho();
+    p_rgh = p - rho*gh;
+}
+else if (pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
 
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-
-Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -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
@@ -29,18 +29,20 @@ Group
 
 Description
     Transient solver for compressible, turbulent flow with a reacting,
-    multiphase particle cloud, and optional sources/constraints.
+    multiphase particle cloud, and surface film modelling.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "basicReactingMultiphaseCloud.H"
+#include "surfaceFilmModel.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "fvOptions.H"
 #include "SLGThermo.H"
+#include "fvOptions.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
@@ -76,10 +78,14 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
 
-        if (!LTS)
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
         {
             #include "compressibleCourantNo.H"
-            #include "setDeltaT.H"
+            #include "setMultiRegionDeltaT.H"
         }
 
         runTime++;
@@ -87,15 +93,16 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         parcels.evolve();
+        surfaceFilm.evolve();
 
-        if (LTS)
+        if (solvePrimaryRegion)
         {
-            #include "setRDeltaT.H"
+            if (pimple.nCorrPIMPLE() <= 1)
+            {
+                #include "rhoEqn.H"
             }
 
-        #include "rhoEqn.H"
+            // --- PIMPLE loop
-
-        // --- Pressure-velocity PIMPLE corrector loop
             while (pimple.loop())
             {
                 #include "UEqn.H"
@@ -115,6 +122,7 @@ int main(int argc, char *argv[])
             }
 
             rho = thermo.rho();
+        }
 
         runTime.write();
 
@@ -123,7 +131,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End\n" << endl;
+    Info<< "End" << endl;
 
     return 0;
 }

applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,15 +36,13 @@ Description
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
+      + surfaceFilm.Srho()
       + fvOptions(rho)
     );
 
     rhoEqn.solve();
 
     fvOptions.correct(rho);
-
-    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
-        << endl;
 }
 
 // ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.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
@@ -33,25 +33,21 @@ Description
 
 if (adjustTimeStep)
 {
-    if (CoNum == -GREAT)
+    const scalar maxDeltaTFact =
-    {
+        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
-        CoNum = SMALL;
+    const scalar deltaTFact =
-    }
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
-
-    const scalar TFactorFluid = maxCo/(CoNum + SMALL);
-    const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
-
-    const scalar dt0 = runTime.deltaTValue();
 
     runTime.setDeltaT
     (
         min
         (
-            dt0*min(min(TFactorFluid, TFactorFilm), 1.2),
+            deltaTFact*runTime.deltaTValue(),
             maxDeltaT
         )
     );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }
 
-
 // ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,9 +24,6 @@ License
 Application
     simpleReactingParcelFoam
 
-Group
-    grpLagrangianSolvers
-
 Description
     Steady state solver for compressible, turbulent flow with reacting,
     multiphase particle clouds and optional sources/constraints.
@@ -38,6 +35,7 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
+#include "IOporosityModelList.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "simpleControl.H"

applications/solvers/lagrangian/sprayFoam/EEqn.H

@@ -19,7 +19,6 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
-      + surfaceFilm.Sh()
       + radiation->Sh(thermo, he)
       + Qdot
       + fvOptions(rho, he)
@@ -36,6 +35,6 @@
     thermo.correct();
     radiation->correct();
 
-    Info<< "T gas min/max   = " << min(T).value() << ", "
+    Info<< "T gas min/max   " << min(T).value() << ", "
         << max(T).value() << endl;
 }

applications/solvers/lagrangian/sprayFoam/YEqn.H

@@ -9,7 +9,6 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 
-
 {
     combustion->correct();
     Qdot = combustion->Qdot();
@@ -25,12 +24,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->alphaEff(), Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
                 parcels.SYi(i, Yi)
-              + fvOptions(rho, Yi)
               + combustion->R(Yi)
-              + surfaceFilm.Srho(i)
+              + fvOptions(rho, Yi)
             );
 
             YEqn.relax();

applications/solvers/lagrangian/sprayFoam/createFieldRefs.H

@@ -1,5 +1,3 @@
-const label inertIndex(composition.species()[inertSpecie]);
-
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-filmModelType& surfaceFilm = tsurfaceFilm();
+const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,9 @@ int main(int argc, char *argv[])
             // Store momentum to set rhoUf for introduced faces.
             volVectorField rhoU("rhoU", rho*U);
 
+            // Store the particle positions
+            parcels.storeGlobalPositions();
+
             // Do any mesh changes
             mesh.update();
 

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files

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

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options

@@ -1,46 +1,36 @@
 EXE_INC = \
-    -I. \
+    -I.. \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I${LIB_SRC}/sampling/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I${LIB_SRC}/meshTools/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
-    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 EXE_LIBS = \
-    -lfiniteVolume \
+    -llagrangian \
-    -lfvOptions \
-    -lsampling \
-    -lmeshTools \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lspecie \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lreactionThermophysicalModels \
-    -lSLGThermo \
-    -lchemistryModel \
-    -lregionModels \
-    -lradiationModels \
-    -lsurfaceFilmModels \
-    -lsurfaceFilmDerivedFvPatchFields \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lODE \
+    -lcompressibleTransportModels \
-    -lcombustionModels
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lradiationModels \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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
+    uncoupledKinematicParcelDyMFoam
+
+Description
+    Transient solver for the passive transport of a particle cloud.
+
+    Uses a pre- calculated velocity field to evolve the cloud.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "basicKinematicCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "cloudName",
+        "name",
+        "specify alternative cloud name. default is 'kinematicCloud'"
+    );
+
+    #define NO_CONTROL
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createFields.H"
+    #include "compressibleCourantNo.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        kinematicCloud.storeGlobalPositions();
+
+        mesh.update();
+
+        U.correctBoundaryConditions();
+
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.evolve();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -118,27 +118,27 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
+        alphaPhi10 = talphaPhi1UD();
 
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
             MULES::correct
             (
                 alphac,
                 alpha1,
-                alphaPhi,
+                alphaPhi10,
-                talphaPhiCorr0.ref(),
+                talphaPhi1Corr0.ref(),
                 zeroField(), zeroField(),
                 1, 0
             );
 
-            alphaPhi += talphaPhiCorr0();
+            alphaPhi10 += talphaPhi1Corr0();
         }
 
         // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
+        talphaPhi1Corr0 = talphaPhi1UD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -152,7 +152,7 @@
 
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
+        tmp<surfaceScalarField> talphaPhi1Un
         (
             fvc::flux
             (
@@ -170,15 +170,15 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 alphac,
                 alpha1,
-                talphaPhiUn(),
+                talphaPhi1Un(),
-                talphaPhiCorr.ref(),
+                talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
                 1,
@@ -188,24 +188,24 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                alphaPhi += talphaPhiCorr();
+                alphaPhi10 += talphaPhi1Corr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
+                alphaPhi10 += 0.5*talphaPhi1Corr();
             }
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi10 = talphaPhi1Un;
 
             MULES::explicitSolve
             (
                 alphac,
                 alpha1,
                 phiCN,
-                alphaPhi,
+                alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
                 1,
@@ -220,12 +220,12 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
-        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
     }
     else
     {
-        talphaPhiCorr0.clear();
+        talphaPhi1Corr0.clear();
     }
 
     if
@@ -235,19 +235,20 @@
     )
     {
         #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
         #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + alphaPhic*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + alphaPhic*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -65,6 +65,14 @@
         phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
     }
 
+    // Add the optional shear compression contribution
+    if (scAlpha > 0)
+    {
+        phic +=
+            scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(U))));
+    }
+
+
     surfaceScalarField::Boundary& phicBf =
         phic.boundaryFieldRef();
 
@@ -105,6 +113,8 @@
                 phiCN,
                 upwind<scalar>(mesh, phiCN)
             ).fvmDiv(phiCN, alpha1)
+       // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
+       //           + fvc::div(phiCN), alpha1)
          ==
             Su + fvm::Sp(Sp + divU, alpha1)
         );
@@ -117,19 +127,19 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
+        alphaPhi10 = talphaPhi1UD();
 
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
+            MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
 
-            alphaPhi += talphaPhiCorr0();
+            alphaPhi10 += talphaPhi1Corr0();
         }
 
         // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
+        talphaPhi1Corr0 = talphaPhi1UD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -143,7 +153,7 @@
 
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
+        tmp<surfaceScalarField> talphaPhi1Un
         (
             fvc::flux
             (
@@ -161,15 +171,15 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 geometricOneField(),
                 alpha1,
-                talphaPhiUn(),
+                talphaPhi1Un(),
-                talphaPhiCorr.ref(),
+                talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
                 1,
@@ -179,24 +189,24 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                alphaPhi += talphaPhiCorr();
+                alphaPhi10 += talphaPhi1Corr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
+                alphaPhi10 += 0.5*talphaPhi1Corr();
             }
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi10 = talphaPhi1Un;
 
             MULES::explicitSolve
             (
                 geometricOneField(),
                 alpha1,
                 phiCN,
-                alphaPhi,
+                alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
                 1,
@@ -211,34 +221,37 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
-        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
     }
     else
     {
-        talphaPhiCorr0.clear();
+        talphaPhi1Corr0.clear();
     }
 
+    #include "rhofs.H"
+
     if
     (
         word(mesh.ddtScheme("ddt(rho,U)"))
      == fv::EulerDdtScheme<vector>::typeName
+     || word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::localEulerDdtScheme<vector>::typeName
     )
     {
-        #include "rhofs.H"
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
-        #include "rhofs.H"
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/VoF/createAlphaFluxes.H

@@ -1,20 +1,21 @@
-IOobject alphaPhiHeader
+IOobject alphaPhi10Header
 (
-    "alphaPhi",
+    "alphaPhi10",
     runTime.timeName(),
     mesh,
     IOobject::READ_IF_PRESENT,
     IOobject::AUTO_WRITE
 );
 
-const bool alphaRestart = alphaPhiHeader.typeHeaderOk<surfaceScalarField>(true);
+const bool alphaRestart =
+    alphaPhi10Header.typeHeaderOk<surfaceScalarField>(true);
 
 // MULES flux from previous time-step
-surfaceScalarField alphaPhi
+surfaceScalarField alphaPhi10
 (
-    alphaPhiHeader,
+    alphaPhi10Header,
     phi*fvc::interpolate(alpha1)
 );
 
 // MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
+tmp<surfaceScalarField> talphaPhi1Corr0;