diff --git a/applications/solvers/compressible/rhoPimpleFoam/Allwmake b/applications/solvers/compressible/rhoPimpleFoam/Allwmake index 241e22eb1b..f6f8ad3635 100755 --- a/applications/solvers/compressible/rhoPimpleFoam/Allwmake +++ b/applications/solvers/compressible/rhoPimpleFoam/Allwmake @@ -3,6 +3,7 @@ cd ${0%/*} || exit 1 # run from this directory set -x wmake +wmake rhoPimplecFoam wmake rhoPorousMRFPimpleFoam wmake rhoPorousMRFLTSPimpleFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H index 1e39e983e4..b7b4b04db5 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H @@ -4,13 +4,12 @@ tmp UEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) + - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U) + turbulence->divDevRhoReff(U) ); UEqn().relax(); -volScalarField rAU(1.0/UEqn().A()); - if (pimple.momentumPredictor()) { solve(UEqn() == -fvc::grad(p)); diff --git a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H index b1d9e4e8b0..57d72887e6 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H @@ -3,10 +3,14 @@ ( fvm::ddt(rho, h) + fvm::div(phi, h) + - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h) - fvm::laplacian(turbulence->alphaEff(), h) == dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) + - ( + fvc::ddt(rho, K) + fvc::div(phi, K) + - (fvc::ddt(rho) + fvc::div(phi))*K + ) ); hEqn.relax(); diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index cfe34f2873..3d324e085b 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -3,6 +3,7 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); +volScalarField rAU(1.0/UEqn().A()); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn().H(); @@ -23,13 +24,15 @@ if (pimple.transonic()) ) ); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); @@ -52,6 +55,8 @@ else ) ); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { // Pressure corrector @@ -59,7 +64,7 @@ else ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C index 26d11a14aa..28787716e0 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,7 +66,10 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files new file mode 100644 index 0000000000..5bb1f8d023 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files @@ -0,0 +1,3 @@ +rhoPimplecFoam.C + +EXE = $(FOAM_APPBIN)/rhoPimplecFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options new file mode 100644 index 0000000000..d16a0ee1d8 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options @@ -0,0 +1,15 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lbasicThermophysicalModels \ + -lspecie \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H new file mode 100644 index 0000000000..7ef1b62533 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H @@ -0,0 +1,117 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn().A()); +volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1())); + +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); + +if (pimple.nCorrPIMPLE() <= 1) +{ + UEqn.clear(); +} + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + surfaceScalarField phic + ( + "phic", + fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() + ); + + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + + fvc::div(phic) + - fvm::laplacian(Dp, p) + ); + + // Relax the pressure equation to maintain diagonal dominance + pEqn.relax(); + + pEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi == phic + pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(Dp, p) + ); + + pEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +U = HbyA - rAtU*fvc::grad(p); +U.correctBoundaryConditions(); +K = 0.5*magSqr(U); + +dpdt = fvc::ddt(p); + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); + +if (!pimple.transonic()) +{ + rho.relax(); +} + +Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C new file mode 100644 index 0000000000..0b01f51c82 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 3 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, see . + +Application + rhoPimplecFoam + +Description + Transient solver for laminar or turbulent flow of compressible fluids + for HVAC and similar applications. + + Uses the flexible PIMPLEC (PISOC-SIMPLEC) solution for time-resolved and + pseudo-transient simulations. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "basicPsiThermo.H" +#include "turbulenceModel.H" +#include "bound.H" +#include "pimpleControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + pimpleControl pimple(mesh); + + #include "createFields.H" + #include "initContinuityErrs.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "hEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C index de32f1f5c6..e305916f29 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,10 @@ int main(int argc, char *argv[]) #include "setrDeltaT.H" - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H index c63f0f8e4b..c24b1f587a 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H @@ -13,8 +13,6 @@ UEqn().relax(); mrfZones.addCoriolis(rho, UEqn()); pZones.addResistance(UEqn()); -volScalarField rAU(1.0/UEqn().A()); - if (pimple.momentumPredictor()) { solve(UEqn() == -fvc::grad(p)); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H index c5c7602a43..3c75e87e64 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H @@ -25,13 +25,15 @@ if (pimple.transonic()) ); mrfZones.relativeFlux(fvc::interpolate(psi), phid); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); @@ -56,6 +58,8 @@ else mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { // Pressure corrector @@ -63,7 +67,7 @@ else ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C index 068de89952..e02787c35d 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,7 +69,10 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H index 27bfb83133..8a7d636f99 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H @@ -20,7 +20,7 @@ ); volScalarField& p = thermo.p(); - volScalarField& h = thermo.h(); + volScalarField& e = thermo.e(); const volScalarField& psi = thermo.psi(); Info<< "Reading field U\n" << endl; diff --git a/applications/solvers/compressible/rhoSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/eEqn.H new file mode 100644 index 0000000000..a1ea771573 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/eEqn.H @@ -0,0 +1,18 @@ +{ + // Kinetic + pressure energy + volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho); + + fvScalarMatrix eEqn + ( + fvm::div(phi, e) + - fvm::Sp(fvc::div(phi), e) + - fvm::laplacian(turbulence->alphaEff(), e) + == + fvc::div(phi)*Ekp - fvc::div(phi, Ekp) + ); + + eEqn.relax(); + eEqn.solve(); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H deleted file mode 100644 index 879578cf49..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H +++ /dev/null @@ -1,18 +0,0 @@ -{ - volScalarField K("K", 0.5*magSqr(U)); - - fvScalarMatrix hEqn - ( - fvm::div(phi, h) - - fvm::Sp(fvc::div(phi), h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - fvc::div(phi)*K - fvc::div(phi, K) - ); - - hEqn.relax(); - - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index cefee48969..6c64056be2 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -54,7 +54,8 @@ else { fvScalarMatrix pEqn ( - fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA) + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) ); pEqn.setReference(pRefCell, pRefValue); @@ -63,7 +64,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi = phiHbyA - pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -88,5 +89,10 @@ if (closedVolume) rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); -rho.relax(); + +if (!simple.transonic()) +{ + rho.relax(); +} + Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H new file mode 100644 index 0000000000..4791062d35 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H @@ -0,0 +1,19 @@ +{ + volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho); + + fvScalarMatrix eEqn + ( + fvm::div(phi, e) + - fvm::Sp(fvc::div(phi), e) + - fvm::laplacian(turbulence->alphaEff(), e) + == + fvc::div(phi)*Ekp - fvc::div(phi, Ekp) + ); + + //pZones.addEnergySource(thermo, rho, eEqn); + + eEqn.relax(); + eEqn.solve(); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/hEqn.H deleted file mode 100644 index cb989b804f..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/hEqn.H +++ /dev/null @@ -1,18 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::div(phi, h) - - fvm::Sp(fvc::div(phi), h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") - ); - - pZones.addEnthalpySource(thermo, rho, hEqn); - - hEqn.relax(); - - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C index 1cb421e5f4..c866aea7df 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -63,7 +63,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "hEqn.H" + #include "eEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C index e87a35fab2..bb249f5430 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,7 +59,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "hEqn.H" + #include "eEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H deleted file mode 100644 index 46a382864e..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H +++ /dev/null @@ -1,61 +0,0 @@ - Info<< "Reading thermophysical properties\n" << endl; - - autoPtr pThermo - ( - basicPsiThermo::New(mesh) - ); - basicPsiThermo& thermo = pThermo(); - - volScalarField rho - ( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - thermo.rho() - ); - - volScalarField& p = thermo.p(); - volScalarField& h = thermo.h(); - const volScalarField& psi = thermo.psi(); - - Info<< "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - #include "compressibleCreatePhi.H" - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, simple.dict(), pRefCell, pRefValue); - - dimensionedScalar rhoMax(simple.dict().lookup("rhoMax")); - dimensionedScalar rhoMin(simple.dict().lookup("rhoMin")); - - Info<< "Creating turbulence model\n" << endl; - autoPtr turbulence - ( - compressible::RASModel::New - ( - rho, - U, - phi, - thermo - ) - ); - - dimensionedScalar initialMass = fvc::domainIntegrate(rho); diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H index 379e720c04..8c7405c346 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H @@ -3,13 +3,11 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); -volScalarField p0(p); - -volScalarField AU(UEqn().A()); -volScalarField AtU(AU - UEqn().H1()); +volScalarField rAU(1.0/UEqn().A()); +volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1())); volVectorField HbyA("HbyA", U); -HbyA = UEqn().H()/AU; +HbyA = rAU*UEqn().H(); UEqn.clear(); @@ -17,28 +15,32 @@ bool closedVolume = false; if (simple.transonic()) { + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + surfaceScalarField phic + ( + "phic", + fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() + ); + + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + while (simple.correctNonOrthogonal()) { - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi*HbyA) & mesh.Sf() - ); - - surfaceScalarField phic - ( - "phic", - fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf() - ); - fvScalarMatrix pEqn ( fvm::div(phid, p) + fvc::div(phic) - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); - // Relax the pressure equation to ensure diagonal-dominance + // Relax the pressure equation to maintain diagonal dominance pEqn.relax(); pEqn.setReference(pRefCell, pRefValue); @@ -53,21 +55,25 @@ if (simple.transonic()) } else { + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + closedVolume = adjustPhi(phiHbyA, U, p); + + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + while (simple.correctNonOrthogonal()) { - surfaceScalarField phiHbyA - ( - "phiHbyA", - fvc::interpolate(rho*HbyA) & mesh.Sf() - ); - - closedVolume = adjustPhi(phi, U, p); - phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf(); - fvScalarMatrix pEqn ( fvc::div(phiHbyA) - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); pEqn.setReference(pRefCell, pRefValue); @@ -81,16 +87,14 @@ else } } -// The incompressibe for of the continuity error check is appropriate for +// The incompressibe form of the continuity error check is appropriate for // steady-state compressible also. #include "incompressible/continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); -U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); -//U = HbyA - fvc::grad(p)/AU; - +U = HbyA - rAtU*fvc::grad(p); U.correctBoundaryConditions(); // For closed-volume cases adjust the pressure and density levels @@ -101,6 +105,7 @@ if (closedVolume) /fvc::domainIntegrate(psi); } +// Recalculate density from the relaxed pressure rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C index 37334e7955..198b5eea40 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C @@ -61,8 +61,8 @@ int main(int argc, char *argv[]) // Velocity-pressure-enthalpy SIMPLEC corrector { #include "UEqn.H" + #include "eEqn.H" #include "pEqn.H" - #include "hEqn.H" } turbulence->correct(); diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index 46cc36216c..aeebee2df5 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -14,13 +14,16 @@ surfaceScalarField phid ) ); + +volScalarField Dp("Dp", rho*rAU); + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H index cad5bb7217..21dc48614e 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H @@ -14,13 +14,15 @@ surfaceScalarField phid ) ); +volScalarField Dp("Dp", rho*rAU); + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C index 6388caea69..19a530a6aa 100644 --- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,13 +85,14 @@ int main(int argc, char *argv[]) ); phi = (rhoO/psi)*phid; + volScalarField Dp("Dp", rho*rAU); fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phi) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H index 102ba89a09..bbbd079ef2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H @@ -9,8 +9,8 @@ tmp tcp = thermo.Cp(); const volScalarField& cp = tcp(); - tmp tkappa = thermo.K(); - //tmp tkappa = thermo.directionalK()*betav; + tmp tkappa = thermo.kappa(); + //tmp tkappa = thermo.directionalKappa()*betav; const volScalarField& kappa = tkappa(); //const volSymmTensorField& K = tK(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H index 7ff98e6211..73269420c0 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H @@ -7,7 +7,7 @@ tmp tcp = thermo.Cp(); const volScalarField& cp = tcp(); - tmp tkappa = thermo.K(); + tmp tkappa = thermo.kappa(); //tmp tkappa = thermo.directionalkappa(); const volScalarField& kappa = tkappa(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H index f4c7081e9b..6f4834000d 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H @@ -9,8 +9,8 @@ tmp tcp = thermo.Cp(); const volScalarField& cp = tcp(); - tmp tkappa = thermo.K(); - //tmp tkappa = thermo.directionalK()*betav; + tmp tkappa = thermo.kappa(); + //tmp tkappa = thermo.directionalKappa()*betav; const volScalarField& kappa = tkappa(); //const volSymmTensorField& K = tK(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H index 8cf6f1c38c..47bd9dad9a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H @@ -7,9 +7,9 @@ tmp tcp = thermo.Cp(); const volScalarField& cp = tcp(); - tmp tkappa = thermo.K(); + tmp tkappa = thermo.kappa(); const volScalarField& kappa = tkappa(); - //tmp tkappa = thermo.directionalK(); + //tmp tkappa = thermo.directionalKappa(); //const volSymmTensorField& kappa = tkappa(); volScalarField& T = thermo.T(); diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files index d89ca6e737..a15838abe8 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/Make/files +++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files @@ -1,3 +1,4 @@ collapseEdges.C +pointEdgeCollapse/pointEdgeCollapse.C EXE = $(FOAM_APPBIN)/collapseEdges diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/options b/applications/utilities/mesh/advanced/collapseEdges/Make/options index 21b17b14c9..d1efa61fd5 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/Make/options +++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options @@ -1,6 +1,10 @@ EXE_INC = \ - -I$(LIB_SRC)/dynamicMesh/lnInclude + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -IpointEdgeCollapse EXE_LIBS = \ -ldynamicMesh \ - -lmeshTools + -lmeshTools \ + -lfiniteVolume diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index f1eb28e594..899159b2a2 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,20 +41,435 @@ Description #include "argList.H" #include "Time.H" -#include "edgeCollapser.H" #include "polyTopoChange.H" #include "polyTopoChanger.H" #include "polyMesh.H" +#include "fvMesh.H" #include "mapPolyMesh.H" #include "mathematicalConstants.H" #include "PackedBoolList.H" #include "SortableList.H" #include "unitConversion.H" +#include "globalMeshData.H" +#include "globalIndex.H" +#include "PointEdgeWave.H" +#include "pointEdgeCollapse.H" +#include "motionSmoother.H" + +#include "OFstream.H" +#include "meshTools.H" using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +label findIndex +( + const labelList& elems, + const label start, + const label size, + const label val +) +{ + for (label i = start; i < size; i++) + { + if (elems[i] == val) + { + return i; + } + } + return -1; +} + + +void filterFace +( + const label faceI, + face& f, + const List& allPointInfo, + const Map >& collapseStrings +) +{ + label newFp = 0; + + face oldFace = f; + + forAll(f, fp) + { + label pointI = f[fp]; + + label collapsePoint = allPointInfo[pointI].collapseIndex(); + + if (collapseStrings.found(collapsePoint)) + { + collapsePoint = collapseStrings[collapsePoint][0]; + } + + if (collapsePoint == -1) + { + WarningIn + ( + "filterFace" + "(const label, face&, const List&)" + ) + << "Point " << pointI << " was not visited by PointEdgeWave" + << endl; + } + else if (collapsePoint == -2) + { + f[newFp++] = pointI; + } + else + { + if (findIndex(f, 0, newFp, collapsePoint) == -1) + { + f[newFp++] = collapsePoint; + } + } + } + + + // Check for pinched face. Tries to correct + // - consecutive duplicate vertex. Removes duplicate vertex. + // - duplicate vertex with one other vertex in between (spike). + // Both of these should not really occur! and should be checked before + // collapsing edges. + + const label size = newFp; + + newFp = 2; + + for (label fp = 2; fp < size; fp++) + { + label fp1 = fp-1; + label fp2 = fp-2; + + label pointI = f[fp]; + + // Search for previous occurrence. + label index = findIndex(f, 0, fp, pointI); + + if (index == fp1) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Removing consecutive duplicate vertex in face " + << f << endl; + // Don't store current pointI + } + else if (index == fp2) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Removing non-consecutive duplicate vertex in face " + << f << endl; + // Don't store current pointI and remove previous + newFp--; + } + else if (index != -1) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Pinched face " << f << endl; + f[newFp++] = pointI; + } + else + { + f[newFp++] = pointI; + } + } + + f.setSize(newFp); +} + + +bool setRefinement +( + const polyMesh& mesh, + polyTopoChange& meshMod, + const List& allPointInfo +) +{ + const cellList& cells = mesh.cells(); + const labelList& faceOwner = mesh.faceOwner(); + const labelList& faceNeighbour = mesh.faceNeighbour(); + const labelListList& pointFaces = mesh.pointFaces(); + const pointZoneMesh& pointZones = mesh.pointZones(); + + globalIndex globalStrings(mesh.nPoints()); + + boolList removedPoints(mesh.nPoints(), false); + + // Create strings of edges + Map > collapseStrings; + + forAll(allPointInfo, pointI) + { + const label collapseIndex = allPointInfo[pointI].collapseIndex(); + + if (collapseIndex != -1 && collapseIndex != -2) + { + collapseStrings(collapseIndex).append(pointI); + } + } + + bool meshChanged = false; + + // Current faces (is also collapseStatus: f.size() < 3) + faceList newFaces(mesh.faces()); + + // Current cellCollapse status + boolList cellRemoved(mesh.nCells(), false); + + label nUnvisited = 0; + label nUncollapsed = 0; + label nCollapsed = 0; + + forAll(allPointInfo, pI) + { + const pointEdgeCollapse& pec = allPointInfo[pI]; + + if (pec.collapseIndex() == -1) + { + nUnvisited++; + } + else if (pec.collapseIndex() == -2) + { + nUncollapsed++; + } + else if (pec.collapseIndex() >= 0) + { + nCollapsed++; + } + } + + label nPoints = allPointInfo.size(); + + reduce(nPoints, sumOp