Transient coupled solvers for incompressible flows

applications/solvers/coupled/transientDyMFoam/CMakeLists.txt

@@ -0,0 +1,46 @@
+# --------------------------------------------------------------------------
+#   ========                 |
+#   \      /  F ield         | foam-extend: Open Source CFD
+#    \    /   O peration     | Version:     4.1
+#     \  /    A nd           | Web:         http://www.foam-extend.org
+#      \/     M anipulation  | For copyright notice see file Copyright
+# --------------------------------------------------------------------------
+# License
+#     This file is part of foam-extend.
+#
+#     foam-extend 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.
+#
+#     foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Description
+#     CMakeLists.txt file for libraries and applications
+#
+# Author
+#     Henrik Rusche, Wikki GmbH, 2017. All rights reserved
+#
+#
+# --------------------------------------------------------------------------
+
+list(APPEND SOURCES
+  MRFPorousFoam.C
+)
+
+# Set minimal environment for external compilation
+if(NOT FOAM_FOUND)
+  cmake_minimum_required(VERSION 2.8)
+  find_package(FOAM REQUIRED)
+endif()
+
+add_foam_executable(MRFPorousFoam
+  DEPENDS incompressibleRASModels
+  SOURCES ${SOURCES}
+)

applications/solvers/coupled/transientDyMFoam/Make/files

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

applications/solvers/coupled/transientDyMFoam/Make/options

@@ -0,0 +1,20 @@
+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/coupled/transientDyMFoam/UEqn.H

@@ -0,0 +1,19 @@
+{
+    // Momentum equation
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U)
+      + fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+
+    rAU = 1.0/UEqn.A();
+
+    // Under-relax momentum. Note this will destroy the H and A
+    UEqn.relax();
+
+    // Insert momentum equation
+    UpEqn.insertEquation(0, UEqn);
+
+#   include "addBlockCoupledBC.H"
+}

applications/solvers/coupled/transientDyMFoam/addBlockCoupledBC.H

@@ -0,0 +1,56 @@
+    // Hack block-coupled boundary conditions: due for rewrite
+    const volScalarField nuEff = turbulence->nuEff();
+
+    forAll (U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].blockCoupled())
+        {
+            // Insert correcting fully implicit coupling coefficient
+
+            const labelList fc = mesh.boundary()[patchI].faceCells();
+
+            const fvPatchVectorField& Up = U.boundaryField()[patchI];
+
+            // Warning: hacked for nuEff in viscosity
+            const scalarField nutpMagSf =
+                nuEff.boundaryField()[patchI]*
+                mesh.magSf().boundaryField()[patchI];
+
+            // Get boundary condition contribution to matrix diagonal
+            tensorField patchDiag =
+                -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf;
+
+            // Get matrix diagonal
+            CoeffField<vector4>::squareTypeField& blockDiag =
+                UpEqn.diag().asSquare();
+
+            forAll (fc, faceI)
+            {
+                blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx();
+                blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy();
+                blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz();
+
+                blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx();
+                blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy();
+                blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz();
+
+                blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx();
+                blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy();
+                blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz();
+            }
+
+         // Get boundary condition contribution to matrix source
+         vectorField patchSource =
+             -Up.blockGradientBoundaryCoeffs()*nutpMagSf;
+
+         // Get matrix source
+         Field<vector4>& blockSource = UpEqn.source();
+
+            forAll (fc, faceI)
+            {
+                blockSource[fc[faceI]](0) -= patchSource[faceI](0);
+                blockSource[fc[faceI]](1) -= patchSource[faceI](1);
+                blockSource[fc[faceI]](2) -= patchSource[faceI](2);
+            }
+        }
+    }

applications/solvers/coupled/transientDyMFoam/boundPU.H

@@ -0,0 +1,32 @@
+{
+    // Bound the pressure
+    dimensionedScalar p1 = min(p);
+    dimensionedScalar p2 = max(p);
+
+    if (p1 < pMin || p2 > pMax)
+    {
+        Info<< "p: " << p1.value() << " " << p2.value()
+            << ".  Bounding." << endl;
+
+        p.max(pMin);
+        p.min(pMax);
+        p.correctBoundaryConditions();
+    }
+
+    // Bound the velocity
+    volScalarField magU = mag(U);
+    dimensionedScalar U1 = max(magU);
+
+    if (U1 > UMax)
+    {
+        Info<< "U: " << U1.value() << ".  Bounding." << endl;
+
+        volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU)
+            + neg(magU - UMax);
+        Ulimiter.max(scalar(0));
+        Ulimiter.min(scalar(1));
+
+        U *= Ulimiter;
+        U.correctBoundaryConditions();
+    }
+}

applications/solvers/coupled/transientDyMFoam/convergenceCheck.H

@@ -0,0 +1,9 @@
+// Check convergence
+if (maxResidual < convergenceCriterion)
+{
+    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
+
+    // Move to the next time-step
+    break;
+}
+

applications/solvers/coupled/transientDyMFoam/correctPhi.H

@@ -0,0 +1,33 @@
+{
+    volScalarField pcorr("pcorr", p);
+    pcorr *= 0;
+
+    // Initialise flux with interpolated velocity
+    phi = fvc::interpolate(U) & mesh.Sf();
+
+    adjustPhi(phi, U, pcorr);
+
+    mesh.schemesDict().setFluxRequired(pcorr.name());
+
+    // while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
+        );
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        // if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+
+        // Fluxes are corrected to absolute velocity and further corrected
+        // later.  HJ, 6/Feb/2009
+    }
+
+#   include "continuityErrs.H"
+}
+

applications/solvers/coupled/transientDyMFoam/couplingTerms.H

@@ -0,0 +1,15 @@
+{
+    // Calculate grad p coupling matrix. Needs to be here if one uses
+    // gradient schemes with limiters.  VV, 9/June/2014
+    BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+    // Calculate div U coupling.  Could be calculated only once since
+    // it is only geometry dependent.  VV, 9/June/2014
+    BlockLduSystem<vector, scalar> UInp(fvm::UDiv(U));
+
+    // Last argument in insertBlockCoupling says if the column direction
+    // should be incremented. This is needed for arbitrary positioning
+    // of U and p in the system. This could be better. VV, 30/April/2014
+    UpEqn.insertBlockCoupling(0, 3, pInU, true);
+    UpEqn.insertBlockCoupling(3, 0, UInp, false);
+}

applications/solvers/coupled/transientDyMFoam/createControls.H

@@ -0,0 +1,10 @@
+#   include "createTimeControls.H"
+
+    Switch correctPhi(false);
+    Switch checkMeshCourantNo(false);
+
+    // Number of outer correctors
+    label nOuterCorrectors = 1;
+
+    label pRefCell = 0;
+    scalar pRefValue = 0;

applications/solvers/coupled/transientDyMFoam/createFields.H

@@ -0,0 +1,70 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info << "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    // Block vector field for velocity (first entry) and pressure (second
+    // entry).
+    Info << "Creating field Up\n" << endl;
+    volVector4Field Up
+    (
+        IOobject
+        (
+            "Up",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector4("zero", dimless, vector4::zero)
+    );
+
+    Info<< "Creating field rAU\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        runTime.deltaT()
+    );
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/coupled/transientDyMFoam/initConvergenceCheck.H

@@ -0,0 +1,5 @@
+    // initialize values for convergence checks
+    BlockSolverPerformance<vector4> residual;
+
+    scalar maxResidual = 0;
+    scalar convergenceCriterion = 0;

applications/solvers/coupled/transientDyMFoam/pEqn.H

@@ -0,0 +1,18 @@
+    // Pressure parts of the continuity equation
+    surfaceScalarField presSource
+    (
+         "presSource",
+         fvc::interpolate(rAU)*
+         (fvc::interpolate(fvc::grad(p)) & mesh.Sf())
+    );
+
+    fvScalarMatrix pEqn
+    (
+      - fvm::laplacian(rAU, p)
+     ==
+      - fvc::div(presSource)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    UpEqn.insertEquation(3, pEqn);

applications/solvers/coupled/transientDyMFoam/readBlockSolverControls.H

@@ -0,0 +1,29 @@
+{
+    dictionary blockSolverDict = mesh.solutionDict().subDict("blockSolver");
+
+    setRefCell
+    (
+        p,
+        blockSolverDict,
+        pRefCell,
+        pRefValue
+    );
+
+    blockSolverDict.readIfPresent
+    (
+        "nOuterCorrectors",
+        nOuterCorrectors
+    );
+
+    blockSolverDict.readIfPresent
+    (
+        "correctPhi",
+        correctPhi
+    );
+
+    blockSolverDict.readIfPresent
+    (
+        "checkMeshCourantNo",
+        checkMeshCourantNo
+    );
+}

applications/solvers/coupled/transientDyMFoam/readFieldBounds.H

@@ -0,0 +1,14 @@
+    // Read field bounds
+    dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
+
+    // Pressure bounds
+    dimensionedScalar pMin("pMin", p.dimensions(), 0);
+    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+
+    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+
+    // Velocity bound
+    dimensionedScalar UMax("UMax", U.dimensions(), 0);
+
+    fieldBounds.lookup(U.name()) >> UMax.value();
+    dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

applications/solvers/coupled/transientDyMFoam/transientDyMFoam.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    transientDyMFoam
+
+Description
+    Transient solver for incompressible, turbulent flow, with implicit
+    coupling between pressure and velocity achieved by fvBlockMatrix.
+
+    Turbulence is solved using the existing turbulence model structure.
+
+    The solver supports dynamic mesh changes
+
+Authors
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fvBlockMatrix.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+#include "dynamicFvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+#   include "initConvergenceCheck.H"
+#   include "createControls.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.run())
+    {
+#       include "readBlockSolverControls.H"
+#       include "readFieldBounds.H"
+#       include "CourantNo.H"
+#       include "setDeltaT.H"
+
+        // Make the fluxes absolute
+        fvc::makeAbsolute(phi, U);
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+        reduce(meshChanged, orOp<bool>());
+
+#       include "volContinuity.H"
+
+        if (correctPhi && meshChanged)
+        {
+            // Fluxes will be corrected to absolute velocity
+            // HJ, 6/Feb/2009
+#           include "correctPhi.H"
+        }
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+        if (mesh.moving() && checkMeshCourantNo)
+        {
+#           include "meshCourantNo.H"
+        }
+
+        if (meshChanged)
+        {
+#           include "CourantNo.H"
+        }
+
+        for (label i = 0; i < nOuterCorrectors; i++)
+        {
+            p.storePrevIter();
+
+            // Initialize the Up block system
+            fvBlockMatrix<vector4> UpEqn(Up);
+
+            // Assemble and insert momentum equation
+#           include "UEqn.H"
+
+            // Assemble and insert pressure equation
+#           include "pEqn.H"
+
+            // Assemble and insert coupling terms
+#           include "couplingTerms.H"
+
+            // Solve the block matrix
+            residual = UpEqn.solve();
+            maxResidual = cmptMax(residual.initialResidual());
+
+            // Retrieve solution
+            UpEqn.retrieveSolution(0, U.internalField());
+            UpEqn.retrieveSolution(3, p.internalField());
+
+            U.correctBoundaryConditions();
+            p.correctBoundaryConditions();
+
+            phi = (fvc::interpolate(U) & mesh.Sf())
+                + pEqn.flux()
+                + presSource;
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phi, U);
+
+#           include "movingMeshContinuityErrs.H"
+
+#           include "boundPU.H"
+
+            p.relax();
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+#       include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}

applications/solvers/coupled/transientFoam/CMakeLists.txt

@@ -0,0 +1,46 @@
+# --------------------------------------------------------------------------
+#   ========                 |
+#   \      /  F ield         | foam-extend: Open Source CFD
+#    \    /   O peration     | Version:     4.1
+#     \  /    A nd           | Web:         http://www.foam-extend.org
+#      \/     M anipulation  | For copyright notice see file Copyright
+# --------------------------------------------------------------------------
+# License
+#     This file is part of foam-extend.
+#
+#     foam-extend 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.
+#
+#     foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Description
+#     CMakeLists.txt file for libraries and applications
+#
+# Author
+#     Henrik Rusche, Wikki GmbH, 2017. All rights reserved
+#
+#
+# --------------------------------------------------------------------------
+
+list(APPEND SOURCES
+  MRFPorousFoam.C
+)
+
+# Set minimal environment for external compilation
+if(NOT FOAM_FOUND)
+  cmake_minimum_required(VERSION 2.8)
+  find_package(FOAM REQUIRED)
+endif()
+
+add_foam_executable(MRFPorousFoam
+  DEPENDS incompressibleRASModels
+  SOURCES ${SOURCES}
+)

applications/solvers/coupled/transientFoam/Make/files

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

applications/solvers/coupled/transientFoam/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/coupled/transientFoam/UEqn.H

@@ -0,0 +1,19 @@
+{
+    // Momentum equation
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U)
+      + fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+
+    rAU = 1.0/UEqn.A();
+
+    // Under-relax momentum. Note this will destroy the H and A
+    UEqn.relax();
+
+    // Insert momentum equation
+    UpEqn.insertEquation(0, UEqn);
+
+#   include "addBlockCoupledBC.H"
+}

applications/solvers/coupled/transientFoam/addBlockCoupledBC.H

@@ -0,0 +1,56 @@
+    // Hack block-coupled boundary conditions: due for rewrite
+    const volScalarField nuEff = turbulence->nuEff();
+
+    forAll (U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].blockCoupled())
+        {
+            // Insert correcting fully implicit coupling coefficient
+
+            const labelList fc = mesh.boundary()[patchI].faceCells();
+
+            const fvPatchVectorField& Up = U.boundaryField()[patchI];
+
+            // Warning: hacked for nuEff in viscosity
+            const scalarField nutpMagSf =
+                nuEff.boundaryField()[patchI]*
+                mesh.magSf().boundaryField()[patchI];
+
+            // Get boundary condition contribution to matrix diagonal
+            tensorField patchDiag =
+                -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf;
+
+            // Get matrix diagonal
+            CoeffField<vector4>::squareTypeField& blockDiag =
+                UpEqn.diag().asSquare();
+
+            forAll (fc, faceI)
+            {
+                blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx();
+                blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy();
+                blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz();
+
+                blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx();
+                blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy();
+                blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz();
+
+                blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx();
+                blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy();
+                blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz();
+            }
+
+         // Get boundary condition contribution to matrix source
+         vectorField patchSource =
+             -Up.blockGradientBoundaryCoeffs()*nutpMagSf;
+
+         // Get matrix source
+         Field<vector4>& blockSource = UpEqn.source();
+
+            forAll (fc, faceI)
+            {
+                blockSource[fc[faceI]](0) -= patchSource[faceI](0);
+                blockSource[fc[faceI]](1) -= patchSource[faceI](1);
+                blockSource[fc[faceI]](2) -= patchSource[faceI](2);
+            }
+        }
+    }

applications/solvers/coupled/transientFoam/boundPU.H

@@ -0,0 +1,32 @@
+{
+    // Bound the pressure
+    dimensionedScalar p1 = min(p);
+    dimensionedScalar p2 = max(p);
+
+    if (p1 < pMin || p2 > pMax)
+    {
+        Info<< "p: " << p1.value() << " " << p2.value()
+            << ".  Bounding." << endl;
+
+        p.max(pMin);
+        p.min(pMax);
+        p.correctBoundaryConditions();
+    }
+
+    // Bound the velocity
+    volScalarField magU = mag(U);
+    dimensionedScalar U1 = max(magU);
+
+    if (U1 > UMax)
+    {
+        Info<< "U: " << U1.value() << ".  Bounding." << endl;
+
+        volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU)
+            + neg(magU - UMax);
+        Ulimiter.max(scalar(0));
+        Ulimiter.min(scalar(1));
+
+        U *= Ulimiter;
+        U.correctBoundaryConditions();
+    }
+}

applications/solvers/coupled/transientFoam/convergenceCheck.H

@@ -0,0 +1,9 @@
+// Check convergence
+if (maxResidual < convergenceCriterion)
+{
+    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
+
+    // Move to the next time-step
+    break;
+}
+

applications/solvers/coupled/transientFoam/couplingTerms.H

@@ -0,0 +1,15 @@
+{
+    // Calculate grad p coupling matrix. Needs to be here if one uses
+    // gradient schemes with limiters.  VV, 9/June/2014
+    BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+    // Calculate div U coupling.  Could be calculated only once since
+    // it is only geometry dependent.  VV, 9/June/2014
+    BlockLduSystem<vector, scalar> UInp(fvm::UDiv(U));
+
+    // Last argument in insertBlockCoupling says if the column direction
+    // should be incremented. This is needed for arbitrary positioning
+    // of U and p in the system. This could be better. VV, 30/April/2014
+    UpEqn.insertBlockCoupling(0, 3, pInU, true);
+    UpEqn.insertBlockCoupling(3, 0, UInp, false);
+}

applications/solvers/coupled/transientFoam/createFields.H

@@ -0,0 +1,70 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info << "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    // Block vector field for velocity (first entry) and pressure (second
+    // entry).
+    Info << "Creating field Up\n" << endl;
+    volVector4Field Up
+    (
+        IOobject
+        (
+            "Up",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector4("zero", dimless, vector4::zero)
+    );
+
+    Info<< "Creating field rAU\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        runTime.deltaT()
+    );
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/coupled/transientFoam/initConvergenceCheck.H

@@ -0,0 +1,5 @@
+    // initialize values for convergence checks
+    BlockSolverPerformance<vector4> residual;
+
+    scalar maxResidual = 0;
+    scalar convergenceCriterion = 0;

applications/solvers/coupled/transientFoam/pEqn.H

@@ -0,0 +1,18 @@
+    // Pressure parts of the continuity equation
+    surfaceScalarField presSource
+    (
+         "presSource",
+         fvc::interpolate(rAU)*
+         (fvc::interpolate(fvc::grad(p)) & mesh.Sf())
+    );
+
+    fvScalarMatrix pEqn
+    (
+      - fvm::laplacian(rAU, p)
+     ==
+      - fvc::div(presSource)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    UpEqn.insertEquation(3, pEqn);

applications/solvers/coupled/transientFoam/readBlockSolverControls.H

@@ -0,0 +1,34 @@
+    label pRefCell = 0;
+    scalar pRefValue = 0;
+
+    // Number of outer correctors
+    const label nOuterCorrectors = readLabel
+    (
+        mesh.solutionDict().subDict("blockSolver").lookup("nOuterCorrectors")
+    );
+
+    setRefCell
+    (
+        p,
+        mesh.solutionDict().subDict("blockSolver"),
+        pRefCell,
+        pRefValue
+    );
+
+    mesh.solutionDict().subDict("blockSolver").readIfPresent
+    (
+        "convergence",
+        convergenceCriterion
+    );
+
+    mesh.solutionDict().subDict("blockSolver").readIfPresent
+    (
+        "convergence",
+        convergenceCriterion
+    );
+
+    mesh.solutionDict().subDict("blockSolver").readIfPresent
+    (
+        "convergence",
+        convergenceCriterion
+    );

applications/solvers/coupled/transientFoam/readFieldBounds.H

@@ -0,0 +1,14 @@
+    // Read field bounds
+    dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
+
+    // Pressure bounds
+    dimensionedScalar pMin("pMin", p.dimensions(), 0);
+    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+
+    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+
+    // Velocity bound
+    dimensionedScalar UMax("UMax", U.dimensions(), 0);
+
+    fieldBounds.lookup(U.name()) >> UMax.value();
+    dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

applications/solvers/coupled/transientFoam/transientFoam.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    transientFoam
+
+Description
+    Transient solver for incompressible, turbulent flow, with implicit
+    coupling between pressure and velocity achieved by fvBlockMatrix.
+    Turbulence is in this version solved using the existing turbulence
+    structure.
+
+Authors
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fvBlockMatrix.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+#   include "initConvergenceCheck.H"
+#   include "createTimeControls.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.run())
+    {
+#       include "readBlockSolverControls.H"
+#       include "readFieldBounds.H"
+
+#       include "CourantNo.H"
+#       include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        for (label i = 0; i < nOuterCorrectors; i++)
+        {
+            p.storePrevIter();
+
+            // Initialize the Up block system
+            fvBlockMatrix<vector4> UpEqn(Up);
+
+            // Assemble and insert momentum equation
+#           include "UEqn.H"
+
+            // Assemble and insert pressure equation
+#           include "pEqn.H"
+
+            // Assemble and insert coupling terms
+#           include "couplingTerms.H"
+
+            // Solve the block matrix
+            residual = UpEqn.solve();
+            maxResidual = cmptMax(residual.initialResidual());
+
+            // Retrieve solution
+            UpEqn.retrieveSolution(0, U.internalField());
+            UpEqn.retrieveSolution(3, p.internalField());
+
+            U.correctBoundaryConditions();
+            p.correctBoundaryConditions();
+
+            phi = (fvc::interpolate(U) & mesh.Sf())
+                + pEqn.flux()
+                + presSource;
+
+
+#           include "continuityErrs.H"
+
+#           include "boundPU.H"
+
+            p.relax();
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+#       include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}