From 76b941e7ff59526f1453305da2fb6cb144334e50 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 27 Apr 2015 15:34:02 +0100 Subject: [PATCH] Coupled MRF porous solver --- .../coupled/MRFPorousFoam/MRFPorousFoam.C | 120 ++++++++++++++++++ .../solvers/coupled/MRFPorousFoam/Make/files | 3 + .../coupled/MRFPorousFoam/Make/options | 12 ++ .../solvers/coupled/MRFPorousFoam/UEqn.H | 15 +++ .../solvers/coupled/MRFPorousFoam/boundPU.H | 32 +++++ .../coupled/MRFPorousFoam/convergenceCheck.H | 8 ++ .../coupled/MRFPorousFoam/couplingTerms.H | 15 +++ .../coupled/MRFPorousFoam/createFields.H | 52 ++++++++ .../solvers/coupled/MRFPorousFoam/createMRF.H | 2 + .../coupled/MRFPorousFoam/createPorous.H | 1 + .../MRFPorousFoam/initConvergenceCheck.H | 4 + .../solvers/coupled/MRFPorousFoam/pEqn.H | 23 ++++ .../MRFPorousFoam/readBlockSolverControls.H | 15 +++ .../coupled/MRFPorousFoam/readFieldBounds.H | 14 ++ 14 files changed, 316 insertions(+) create mode 100644 applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C create mode 100644 applications/solvers/coupled/MRFPorousFoam/Make/files create mode 100644 applications/solvers/coupled/MRFPorousFoam/Make/options create mode 100644 applications/solvers/coupled/MRFPorousFoam/UEqn.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/boundPU.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/couplingTerms.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createFields.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createMRF.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/createPorous.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/pEqn.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C new file mode 100644 index 000000000..9ee0e7f41 --- /dev/null +++ b/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 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; +} diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/files b/applications/solvers/coupled/MRFPorousFoam/Make/files new file mode 100644 index 000000000..2ae8a1e4f --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/Make/files @@ -0,0 +1,3 @@ +MRFPorousFoam.C + +EXE = $(FOAM_APPBIN)/MRFPorousFoam diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/options b/applications/solvers/coupled/MRFPorousFoam/Make/options new file mode 100644 index 000000000..5d2b99e6d --- /dev/null +++ b/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 diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H new file mode 100644 index 000000000..ed0e2bc44 --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/MRFPorousFoam/boundPU.H b/applications/solvers/coupled/MRFPorousFoam/boundPU.H new file mode 100644 index 000000000..eca655e8c --- /dev/null +++ b/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(); + } +} diff --git a/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H new file mode 100644 index 000000000..9c7734131 --- /dev/null +++ b/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; +} + diff --git a/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H b/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/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 pInU(fvm::grad(p)); + + // Calculate div U coupling. Could be calculated only once since + // it is only geometry dependent. VV, 9/June/2014 + BlockLduSystem 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); +} diff --git a/applications/solvers/coupled/MRFPorousFoam/createFields.H b/applications/solvers/coupled/MRFPorousFoam/createFields.H new file mode 100644 index 000000000..51403c4a7 --- /dev/null +++ b/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 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) +); diff --git a/applications/solvers/coupled/MRFPorousFoam/createMRF.H b/applications/solvers/coupled/MRFPorousFoam/createMRF.H new file mode 100644 index 000000000..161446a8e --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createMRF.H @@ -0,0 +1,2 @@ + MRFZones mrfZones(mesh); + mrfZones.correctBoundaryVelocity(U); diff --git a/applications/solvers/coupled/MRFPorousFoam/createPorous.H b/applications/solvers/coupled/MRFPorousFoam/createPorous.H new file mode 100644 index 000000000..430b466aa --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createPorous.H @@ -0,0 +1 @@ + porousZones pZones(mesh); diff --git a/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H new file mode 100644 index 000000000..482bc421d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H @@ -0,0 +1,4 @@ +// initialize values for convergence checks + + scalar maxResidual = 0; + scalar convergenceCriterion = 0; diff --git a/applications/solvers/coupled/MRFPorousFoam/pEqn.H b/applications/solvers/coupled/MRFPorousFoam/pEqn.H new file mode 100644 index 000000000..9029b873f --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H new file mode 100644 index 000000000..93c2bda47 --- /dev/null +++ b/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 + ); diff --git a/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H b/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H new file mode 100644 index 000000000..4a6c6f787 --- /dev/null +++ b/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);