diff --git a/applications/solvers/coupled/transientDyMFoam/CMakeLists.txt b/applications/solvers/coupled/transientDyMFoam/CMakeLists.txt new file mode 100644 index 000000000..1593c7ea2 --- /dev/null +++ b/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 . +# +# 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} +) diff --git a/applications/solvers/coupled/transientDyMFoam/Make/files b/applications/solvers/coupled/transientDyMFoam/Make/files new file mode 100644 index 000000000..a5089d091 --- /dev/null +++ b/applications/solvers/coupled/transientDyMFoam/Make/files @@ -0,0 +1,3 @@ +transientDyMFoam.C + +EXE = $(FOAM_APPBIN)/transientDyMFoam diff --git a/applications/solvers/coupled/transientDyMFoam/Make/options b/applications/solvers/coupled/transientDyMFoam/Make/options new file mode 100644 index 000000000..7b3201e92 --- /dev/null +++ b/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 diff --git a/applications/solvers/coupled/transientDyMFoam/UEqn.H b/applications/solvers/coupled/transientDyMFoam/UEqn.H new file mode 100644 index 000000000..68a917237 --- /dev/null +++ b/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" +} diff --git a/applications/solvers/coupled/transientDyMFoam/addBlockCoupledBC.H b/applications/solvers/coupled/transientDyMFoam/addBlockCoupledBC.H new file mode 100644 index 000000000..b80aeb90e --- /dev/null +++ b/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::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& 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); + } + } + } diff --git a/applications/solvers/coupled/transientDyMFoam/boundPU.H b/applications/solvers/coupled/transientDyMFoam/boundPU.H new file mode 100644 index 000000000..eca655e8c --- /dev/null +++ b/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(); + } +} diff --git a/applications/solvers/coupled/transientDyMFoam/convergenceCheck.H b/applications/solvers/coupled/transientDyMFoam/convergenceCheck.H new file mode 100644 index 000000000..e49880fe7 --- /dev/null +++ b/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; +} + diff --git a/applications/solvers/coupled/transientDyMFoam/correctPhi.H b/applications/solvers/coupled/transientDyMFoam/correctPhi.H new file mode 100644 index 000000000..58573813e --- /dev/null +++ b/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" +} + diff --git a/applications/solvers/coupled/transientDyMFoam/couplingTerms.H b/applications/solvers/coupled/transientDyMFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/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 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/transientDyMFoam/createControls.H b/applications/solvers/coupled/transientDyMFoam/createControls.H new file mode 100644 index 000000000..bff565cce --- /dev/null +++ b/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; diff --git a/applications/solvers/coupled/transientDyMFoam/createFields.H b/applications/solvers/coupled/transientDyMFoam/createFields.H new file mode 100644 index 000000000..0c4fe2b44 --- /dev/null +++ b/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 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()); diff --git a/applications/solvers/coupled/transientDyMFoam/initConvergenceCheck.H b/applications/solvers/coupled/transientDyMFoam/initConvergenceCheck.H new file mode 100644 index 000000000..6b4622db6 --- /dev/null +++ b/applications/solvers/coupled/transientDyMFoam/initConvergenceCheck.H @@ -0,0 +1,5 @@ + // initialize values for convergence checks + BlockSolverPerformance residual; + + scalar maxResidual = 0; + scalar convergenceCriterion = 0; diff --git a/applications/solvers/coupled/transientDyMFoam/pEqn.H b/applications/solvers/coupled/transientDyMFoam/pEqn.H new file mode 100644 index 000000000..b38548b8e --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/transientDyMFoam/readBlockSolverControls.H b/applications/solvers/coupled/transientDyMFoam/readBlockSolverControls.H new file mode 100644 index 000000000..145c4061c --- /dev/null +++ b/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 + ); +} diff --git a/applications/solvers/coupled/transientDyMFoam/readFieldBounds.H b/applications/solvers/coupled/transientDyMFoam/readFieldBounds.H new file mode 100644 index 000000000..4a6c6f787 --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/transientDyMFoam/transientDyMFoam.C b/applications/solvers/coupled/transientDyMFoam/transientDyMFoam.C new file mode 100644 index 000000000..9ba963c8a --- /dev/null +++ b/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 . + +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()); + +# 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 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; +} diff --git a/applications/solvers/coupled/transientFoam/CMakeLists.txt b/applications/solvers/coupled/transientFoam/CMakeLists.txt new file mode 100644 index 000000000..1593c7ea2 --- /dev/null +++ b/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 . +# +# 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} +) diff --git a/applications/solvers/coupled/transientFoam/Make/files b/applications/solvers/coupled/transientFoam/Make/files new file mode 100644 index 000000000..af55c51dd --- /dev/null +++ b/applications/solvers/coupled/transientFoam/Make/files @@ -0,0 +1,3 @@ +transientFoam.C + +EXE = $(FOAM_APPBIN)/transientFoam diff --git a/applications/solvers/coupled/transientFoam/Make/options b/applications/solvers/coupled/transientFoam/Make/options new file mode 100644 index 000000000..ee84294bc --- /dev/null +++ b/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 diff --git a/applications/solvers/coupled/transientFoam/UEqn.H b/applications/solvers/coupled/transientFoam/UEqn.H new file mode 100644 index 000000000..68a917237 --- /dev/null +++ b/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" +} diff --git a/applications/solvers/coupled/transientFoam/addBlockCoupledBC.H b/applications/solvers/coupled/transientFoam/addBlockCoupledBC.H new file mode 100644 index 000000000..b80aeb90e --- /dev/null +++ b/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::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& 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); + } + } + } diff --git a/applications/solvers/coupled/transientFoam/boundPU.H b/applications/solvers/coupled/transientFoam/boundPU.H new file mode 100644 index 000000000..eca655e8c --- /dev/null +++ b/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(); + } +} diff --git a/applications/solvers/coupled/transientFoam/convergenceCheck.H b/applications/solvers/coupled/transientFoam/convergenceCheck.H new file mode 100644 index 000000000..e49880fe7 --- /dev/null +++ b/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; +} + diff --git a/applications/solvers/coupled/transientFoam/couplingTerms.H b/applications/solvers/coupled/transientFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/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 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/transientFoam/createFields.H b/applications/solvers/coupled/transientFoam/createFields.H new file mode 100644 index 000000000..0c4fe2b44 --- /dev/null +++ b/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 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()); diff --git a/applications/solvers/coupled/transientFoam/initConvergenceCheck.H b/applications/solvers/coupled/transientFoam/initConvergenceCheck.H new file mode 100644 index 000000000..6b4622db6 --- /dev/null +++ b/applications/solvers/coupled/transientFoam/initConvergenceCheck.H @@ -0,0 +1,5 @@ + // initialize values for convergence checks + BlockSolverPerformance residual; + + scalar maxResidual = 0; + scalar convergenceCriterion = 0; diff --git a/applications/solvers/coupled/transientFoam/pEqn.H b/applications/solvers/coupled/transientFoam/pEqn.H new file mode 100644 index 000000000..b38548b8e --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/transientFoam/readBlockSolverControls.H b/applications/solvers/coupled/transientFoam/readBlockSolverControls.H new file mode 100644 index 000000000..c8cda55ab --- /dev/null +++ b/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 + ); diff --git a/applications/solvers/coupled/transientFoam/readFieldBounds.H b/applications/solvers/coupled/transientFoam/readFieldBounds.H new file mode 100644 index 000000000..4a6c6f787 --- /dev/null +++ b/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); diff --git a/applications/solvers/coupled/transientFoam/transientFoam.C b/applications/solvers/coupled/transientFoam/transientFoam.C new file mode 100644 index 000000000..012e1d57f --- /dev/null +++ b/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 . + +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 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; +}