From 2d5052d7148ad1537930cdcd694d10471ccfa035 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 16 Mar 2016 17:49:11 +0000 Subject: [PATCH] Temporary: consistent solver examples --- .../consistentIcoFoam/Make/files | 3 + .../consistentIcoFoam/Make/options | 6 + .../consistentIcoFoam/consistentIcoFoam.C | 125 +++++++++++++++++ .../consistentIcoFoam/createFields.H | 55 ++++++++ .../consistentPimpleDyMFoam/Make/files | 3 + .../consistentPimpleDyMFoam/Make/options | 21 +++ .../consistentPimpleDyMFoam/UEqn.H | 33 +++++ .../consistentPimpleDyMFoam.C | 128 ++++++++++++++++++ .../consistentPimpleDyMFoam/correctPhi.H | 58 ++++++++ .../consistentPimpleDyMFoam/createFields.H | 58 ++++++++ .../consistentPimpleDyMFoam/pEqn.H | 71 ++++++++++ .../consistentPimpleDyMFoam/readControls.H | 14 ++ .../consistentSimpleFoam/Make/files | 3 + .../consistentSimpleFoam/Make/options | 12 ++ .../consistentSimpleFoam/UEqn.H | 20 +++ .../consistentSimpleFoam.C | 86 ++++++++++++ .../consistentSimpleFoam/convergenceCheck.H | 9 ++ .../consistentSimpleFoam/createFields.H | 42 ++++++ .../initConvergenceCheck.H | 7 + .../consistentSimpleFoam/pEqn.H | 59 ++++++++ 20 files changed, 813 insertions(+) create mode 100644 applications/solvers/incompressible/consistentIcoFoam/Make/files create mode 100644 applications/solvers/incompressible/consistentIcoFoam/Make/options create mode 100644 applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C create mode 100644 applications/solvers/incompressible/consistentIcoFoam/createFields.H create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/Make/files create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/Make/options create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/UEqn.H create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/consistentPimpleDyMFoam.C create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/correctPhi.H create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/createFields.H create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/pEqn.H create mode 100644 applications/solvers/incompressible/consistentPimpleDyMFoam/readControls.H create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/Make/files create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/Make/options create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/UEqn.H create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/consistentSimpleFoam.C create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/convergenceCheck.H create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/createFields.H create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/initConvergenceCheck.H create mode 100644 applications/solvers/incompressible/consistentSimpleFoam/pEqn.H diff --git a/applications/solvers/incompressible/consistentIcoFoam/Make/files b/applications/solvers/incompressible/consistentIcoFoam/Make/files new file mode 100644 index 000000000..2f6e0d6e9 --- /dev/null +++ b/applications/solvers/incompressible/consistentIcoFoam/Make/files @@ -0,0 +1,3 @@ +consistentIcoFoam.C + +EXE = $(FOAM_APPBIN)/consistentIcoFoam diff --git a/applications/solvers/incompressible/consistentIcoFoam/Make/options b/applications/solvers/incompressible/consistentIcoFoam/Make/options new file mode 100644 index 000000000..552db34e2 --- /dev/null +++ b/applications/solvers/incompressible/consistentIcoFoam/Make/options @@ -0,0 +1,6 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -llduSolvers diff --git a/applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C b/applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C new file mode 100644 index 000000000..459adc717 --- /dev/null +++ b/applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 + consistentIcoFoam + +Description + Transient solver for incompressible, laminar flow of Newtonian fluids. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" + +# include "createTime.H" +# include "createMesh.H" +# include "createFields.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + +# include "readPISOControls.H" +# include "CourantNo.H" + + // Convection-diffusion matrix + fvVectorMatrix HUEqn + ( + fvm::div(phi, U) + - fvm::laplacian(nu, U) + ); + + // ddt matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + solve(ddtUEqn + HUEqn == -fvc::grad(p)); + + // Prepare clean Ap without time derivative contribution + // HJ, 26/Oct/2015 + volScalarField aU = HUEqn.A(); + + // --- PISO loop + + for (int corr = 0; corr < nCorr; corr++) + { + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); + + adjustPhi(phi, U, p); + + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + } + +# include "continuityErrs.H" + + // Note: cannot call H(U) here because the velocity is not complete + // HJ, 22/Jan/2016 + U = 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ); + U.correctBoundaryConditions(); + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/incompressible/consistentIcoFoam/createFields.H b/applications/solvers/incompressible/consistentIcoFoam/createFields.H new file mode 100644 index 000000000..6a7d6b9a8 --- /dev/null +++ b/applications/solvers/incompressible/consistentIcoFoam/createFields.H @@ -0,0 +1,55 @@ + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + dimensionedScalar nu + ( + transportProperties.lookup("nu") + ); + + 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" + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/files b/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/files new file mode 100644 index 000000000..8b528f270 --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/files @@ -0,0 +1,3 @@ +consistentPimpleDyMFoam.C + +EXE = $(FOAM_APPBIN)/consistentPimpleDyMFoam diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/options b/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/options new file mode 100644 index 000000000..3f3ed2223 --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/Make/options @@ -0,0 +1,21 @@ +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 \ + -L$(MESQUITE_LIB_DIR) -lmesquite diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/consistentPimpleDyMFoam/UEqn.H new file mode 100644 index 000000000..2a470de9c --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/UEqn.H @@ -0,0 +1,33 @@ + // Convection-diffusion matrix + fvVectorMatrix HUEqn + ( + fvm::div(phi, U) + + turbulence->divDevReff(U) + ); + + // ddt matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + // Get under-relaxation factor and under-relax the diagonal + scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); + + if (oCorr == nOuterCorr - 1) + { + if (mesh.solutionDict().relax("UFinal")) + { + UUrf = mesh.solutionDict().relaxationFactor("UFinal"); + } + else + { + UUrf = 1; + } + } + + // Solve momentum predictor + solve + ( + ddtUEqn + + relax(HUEqn, UUrf) + == + -fvc::grad(p) + ); diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/consistentPimpleDyMFoam.C b/applications/solvers/incompressible/consistentPimpleDyMFoam/consistentPimpleDyMFoam.C new file mode 100644 index 000000000..5e9e13fbd --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/consistentPimpleDyMFoam.C @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 + pimpleDyMFoam.C + +Description + Transient solver for incompressible, flow of Newtonian fluids + with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.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 "readPIMPLEControls.H" +# include "initContinuityErrs.H" +# include "createFields.H" +# include "readTimeControls.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readControls.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(); + +# include "volContinuity.H" + + if (checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + + // Mesh motion update + if (correctPhi && meshChanged) + { + // Fluxes will be corrected to absolute velocity + // HJ, 6/Feb/2009 +# include "correctPhi.H" + } + + if (meshChanged) + { +# include "CourantNo.H" + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + // --- PIMPLE loop + label oCorr = 0; + do + { + if (nOuterCorr != 1) + { + p.storePrevIter(); + } + +# include "UEqn.H" + + // --- PISO loop + for (int corr = 0; corr < nCorr; corr++) + { +# include "pEqn.H" + } + + turbulence->correct(); + } while (++oCorr < nOuterCorr); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/consistentPimpleDyMFoam/correctPhi.H new file mode 100644 index 000000000..60b36ceef --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/correctPhi.H @@ -0,0 +1,58 @@ +{ +# include "continuityErrs.H" + + wordList pcorrTypes + ( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + for (label i = 0; i < p.boundaryField().size(); i++) + { + if (p.boundaryField()[i].fixesValue()) + { + pcorrTypes[i] = fixedValueFvPatchScalarField::typeName; + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0), + pcorrTypes + ); + + // Initialise flux with interpolated velocity + phi = fvc::interpolate(U) & mesh.Sf(); + + adjustPhi(phi, U, pcorr); + + for(int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(1/aU, pcorr) == fvc::div(phi) + ); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + 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/incompressible/consistentPimpleDyMFoam/createFields.H b/applications/solvers/incompressible/consistentPimpleDyMFoam/createFields.H new file mode 100644 index 000000000..1c04f5dbe --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/createFields.H @@ -0,0 +1,58 @@ + 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" + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue); + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + + Info<< "Reading field aU if present\n" << endl; + volScalarField aU + ( + IOobject + ( + "aU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + 1/runTime.deltaT(), + zeroGradientFvPatchScalarField::typeName + ); diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/consistentPimpleDyMFoam/pEqn.H new file mode 100644 index 000000000..c90d2f5cc --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/pEqn.H @@ -0,0 +1,71 @@ +{ + p.boundaryField().updateCoeffs(); + + // Prepare clean Ap without time derivative contribution and + // without contribution from under-relaxation + // HJ, 26/Oct/2015 + aU = HUEqn.A(); + + // Store velocity under-relaxation point before using U for + // the flux precursor + U.storePrevIter(); + + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); + + // Global flux balance + adjustPhi(phi, U, p); + + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + + if + ( +// oCorr == nOuterCorr - 1 + corr == nCorr - 1 + && nonOrth == nNonOrthCorr + ) + { + pEqn.solve + ( + mesh.solutionDict().solver(p.name() + "Final") + ); + } + else + { + pEqn.solve(mesh.solutionDict().solver(p.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + } + + // Explicitly relax pressure for momentum corrector + if (oCorr != nOuterCorr - 1) + { + p.relax(); + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "movingMeshContinuityErrs.H" + + U = UUrf* + ( + 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ) + ) + + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/incompressible/consistentPimpleDyMFoam/readControls.H b/applications/solvers/incompressible/consistentPimpleDyMFoam/readControls.H new file mode 100644 index 000000000..3bd20c5c5 --- /dev/null +++ b/applications/solvers/incompressible/consistentPimpleDyMFoam/readControls.H @@ -0,0 +1,14 @@ +# include "readTimeControls.H" +# include "readPIMPLEControls.H" + + bool correctPhi = false; + if (pimple.found("correctPhi")) + { + correctPhi = Switch(pimple.lookup("correctPhi")); + } + + bool checkMeshCourantNo = false; + if (pimple.found("checkMeshCourantNo")) + { + checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo")); + } diff --git a/applications/solvers/incompressible/consistentSimpleFoam/Make/files b/applications/solvers/incompressible/consistentSimpleFoam/Make/files new file mode 100644 index 000000000..f11b16baf --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/Make/files @@ -0,0 +1,3 @@ +consistentSimpleFoam.C + +EXE = $(FOAM_APPBIN)/consistentSimpleFoam diff --git a/applications/solvers/incompressible/consistentSimpleFoam/Make/options b/applications/solvers/incompressible/consistentSimpleFoam/Make/options new file mode 100644 index 000000000..3563afc4f --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/Make/options @@ -0,0 +1,12 @@ +EXE_INC = -g \ + -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/incompressible/consistentSimpleFoam/UEqn.H b/applications/solvers/incompressible/consistentSimpleFoam/UEqn.H new file mode 100644 index 000000000..fa9e41162 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/UEqn.H @@ -0,0 +1,20 @@ + // Solve the momentum equation + + tmp HUEqn + ( + fvm::div(phi, U) + + turbulence->divDevReff(U) + ); + + // Get under-relaxation factor and under-relax the diagonal + const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); + + // Momentum solution + eqnResidual = solve + ( + relax(HUEqn(), UUrf) + == + -fvc::grad(p) + ).initialResidual(); + + maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/incompressible/consistentSimpleFoam/consistentSimpleFoam.C b/applications/solvers/incompressible/consistentSimpleFoam/consistentSimpleFoam.C new file mode 100644 index 000000000..9e2f11757 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/consistentSimpleFoam.C @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 + consistentSimpleFoam + +Description + Steady-state solver for incompressible, turbulent flow + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "createFields.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + +# include "readSIMPLEControls.H" +# include "initConvergenceCheck.H" + + p.storePrevIter(); + + // Pressure-velocity SIMPLE corrector + { +# include "UEqn.H" +# include "pEqn.H" + } + + 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/incompressible/consistentSimpleFoam/convergenceCheck.H b/applications/solvers/incompressible/consistentSimpleFoam/convergenceCheck.H new file mode 100644 index 000000000..895806319 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/convergenceCheck.H @@ -0,0 +1,9 @@ +// check convergence + +if (maxResidual < convergenceCriterion) +{ + Info<< "reached convergence criterion: " << convergenceCriterion << endl; + runTime.writeAndEnd(); + Info<< "latestTime = " << runTime.timeName() << endl; +} + diff --git a/applications/solvers/incompressible/consistentSimpleFoam/createFields.H b/applications/solvers/incompressible/consistentSimpleFoam/createFields.H new file mode 100644 index 000000000..ab4249625 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/createFields.H @@ -0,0 +1,42 @@ + 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" + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue); + + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::RASModel::New(U, phi, laminarTransport) + ); diff --git a/applications/solvers/incompressible/consistentSimpleFoam/initConvergenceCheck.H b/applications/solvers/incompressible/consistentSimpleFoam/initConvergenceCheck.H new file mode 100644 index 000000000..b56197f22 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/initConvergenceCheck.H @@ -0,0 +1,7 @@ +// initialize values for convergence checks + + scalar eqnResidual = 1, maxResidual = 0; + scalar convergenceCriterion = 0; + + simple.readIfPresent("convergence", convergenceCriterion); + diff --git a/applications/solvers/incompressible/consistentSimpleFoam/pEqn.H b/applications/solvers/incompressible/consistentSimpleFoam/pEqn.H new file mode 100644 index 000000000..021350f09 --- /dev/null +++ b/applications/solvers/incompressible/consistentSimpleFoam/pEqn.H @@ -0,0 +1,59 @@ + p.boundaryField().updateCoeffs(); + + // Prepare clean 1/Ap without contribution from under-relaxation + // HJ, 26/Oct/2015 + volScalarField rUA + ( + "(1|A(U))", + 1/HUEqn().A() + ); + + // Store velocity under-relaxation point before using U for + // the flux precursor + U.storePrevIter(); + + U = rUA*HUEqn().H(); + HUEqn.clear(); + phi = fvc::interpolate(U) & mesh.Sf(); + + // Global flux balance + adjustPhi(phi, U, p); + + // Non-orthogonal pressure corrector loop + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rUA, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + + // Retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + } + +# include "continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + // Momentum corrector + // Note: since under-relaxation does not change aU, H/a in U can be + // re-used. HJ, 22/Jan/2016 + U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +