Coupled MRF porous solver

applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    MRFPorousFoam
+
+Description
+    Steady-state 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.
+
+    Added support for MRF and porous zones
+
+Authors
+    Hrvoje Jasak, Wikki Ltd.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fvBlockMatrix.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "MRFZones.H"
+#include "porousZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "createFields.H"
+#   include "createMRF.H"
+#   include "createPorous.H"
+#   include "initContinuityErrs.H"
+#   include "initConvergenceCheck.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+#       include "readBlockSolverControls.H"
+#       include "readFieldBounds.H"
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        p.storePrevIter();
+
+        // Initialize the Up block system (matrix, source and reference to Up)
+        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
+        maxResidual = cmptMax(UpEqn.solve().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"
+
+        // Make flux relative in rotating zones
+        mrfZones.relativeFlux(phi);
+
+#       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;
+}

applications/solvers/coupled/MRFPorousFoam/Make/files

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

applications/solvers/coupled/MRFPorousFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -0,0 +1,15 @@
+    // Momentum equation
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U)
+      + fvm::div(phi, U)
+      + turbulence->divDevReff(U)
+    );
+
+    // Add MRF and porous sources
+    mrfZones.addCoriolis(UEqn);
+    pZones.addResistance(UEqn);
+
+    UEqn.relax();
+
+    UpEqn.insertEquation(0, UEqn);

applications/solvers/coupled/MRFPorousFoam/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/MRFPorousFoam/convergenceCheck.H

@@ -0,0 +1,8 @@
+// Check convergence
+if (maxResidual < convergenceCriterion)
+{
+    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
+    runTime.writeAndEnd();
+    Info<< "latestTime = " << runTime.timeName() << endl;
+}
+

applications/solvers/coupled/MRFPorousFoam/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/MRFPorousFoam/createFields.H

@@ -0,0 +1,52 @@
+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::RASModel> turbulence
+(
+    incompressible::RASModel::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)
+);

applications/solvers/coupled/MRFPorousFoam/createMRF.H

@@ -0,0 +1,2 @@
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);

applications/solvers/coupled/MRFPorousFoam/createPorous.H

@@ -0,0 +1 @@
+    porousZones pZones(mesh);

applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H

@@ -0,0 +1,4 @@
+// initialize values for convergence checks
+
+    scalar maxResidual = 0;
+    scalar convergenceCriterion = 0;

applications/solvers/coupled/MRFPorousFoam/pEqn.H

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

applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H

@@ -0,0 +1,15 @@
+    label pRefCell = 0;
+    scalar pRefValue = 0;
+    setRefCell
+    (
+        p,
+        mesh.solutionDict().subDict("blockSolver"),
+        pRefCell,
+        pRefValue
+    );
+
+    mesh.solutionDict().subDict("blockSolver").readIfPresent
+    (
+        "convergence",
+        convergenceCriterion
+    );

applications/solvers/coupled/MRFPorousFoam/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);