diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C index 0733802b08..08073398a7 100644 --- a/applications/solvers/DNS/dnsFoam/dnsFoam.C +++ b/applications/solvers/DNS/dnsFoam/dnsFoam.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,7 +85,7 @@ int main(int argc, char *argv[]) for (int corr=1; corr<=1; corr++) { - volScalarField rAU(1.0/UEqn.A()); + volScalarField rAU("Dp", 1.0/UEqn.A()); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); diff --git a/applications/solvers/combustion/chemFoam/createSingleCellMesh.H b/applications/solvers/combustion/chemFoam/createSingleCellMesh.H index b5a9d283ac..c125a7f31f 100644 --- a/applications/solvers/combustion/chemFoam/createSingleCellMesh.H +++ b/applications/solvers/combustion/chemFoam/createSingleCellMesh.H @@ -23,7 +23,7 @@ fvMesh mesh fvMesh::defaultRegion, runTime.timeName(), runTime, - IOobject::MUST_READ + IOobject::READ_IF_PRESENT ), xferMove >(points), faces.xfer(), diff --git a/applications/solvers/combustion/chemFoam/readInitialConditions.H b/applications/solvers/combustion/chemFoam/readInitialConditions.H index 3fbb13a0fc..08050be816 100644 --- a/applications/solvers/combustion/chemFoam/readInitialConditions.H +++ b/applications/solvers/combustion/chemFoam/readInitialConditions.H @@ -99,7 +99,7 @@ scalar u0 = hs0 - p0/rho0; scalar R0 = p0/(rho0*T0); Rspecific[0] = R0; - + scalar integratedHeat = 0.0; Info << constProp << " will be held constant." << nl diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H index f5364ae314..3f4f485cf6 100644 --- a/applications/solvers/combustion/fireFoam/pEqn.H +++ b/applications/solvers/combustion/fireFoam/pEqn.H @@ -1,11 +1,11 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); +surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); -surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); +surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); surfaceScalarField phiHbyA ( @@ -15,7 +15,7 @@ surfaceScalarField phiHbyA (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) - - phig + + phig ); @@ -37,7 +37,7 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA + p_rghEqn.flux(); - U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/compressible/rhoCentralFoam/readThermophysicalProperties.H b/applications/solvers/compressible/rhoCentralFoam/readThermophysicalProperties.H index f9e309ffec..e0b16f0f39 100644 --- a/applications/solvers/compressible/rhoCentralFoam/readThermophysicalProperties.H +++ b/applications/solvers/compressible/rhoCentralFoam/readThermophysicalProperties.H @@ -16,11 +16,8 @@ IOdictionary thermophysicalProperties dimensionedScalar Pr ( - dimensionedScalar::lookupOrDefault - ( - "Pr", - thermophysicalProperties, - 1.0 - ) + "Pr", + dimless, + thermophysicalProperties.subDict("mixture").subDict("transport") + .lookup("Pr") ); - 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/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C similarity index 68% rename from applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C rename to applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C index 357c7e413a..0b01f51c82 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.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 @@ -22,88 +22,74 @@ License along with OpenFOAM. If not, see . Application - twoPhaseEulerFoam + rhoPimplecFoam Description - Solver for a system of 2 incompressible fluid phases with one phase - dispersed, e.g. gas bubbles in a liquid or solid particles in a gas. + 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 "nearWallDist.H" -#include "wallFvPatch.H" -#include "Switch.H" - -#include "IFstream.H" -#include "OFstream.H" - -#include "dragModel.H" -#include "phaseModel.H" -#include "kineticTheoryModel.H" - +#include "basicPsiThermo.H" +#include "turbulenceModel.H" +#include "bound.H" #include "pimpleControl.H" -#include "MRFZones.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" - #include "createTime.H" #include "createMesh.H" - #include "readGravitationalAcceleration.H" - #include "createFields.H" - #include "readPPProperties.H" - #include "initContinuityErrs.H" - #include "createMRFZones.H" - #include "readTimeControls.H" - #include "CourantNo.H" - #include "setInitialDeltaT.H" pimpleControl pimple(mesh); + #include "createFields.H" + #include "initContinuityErrs.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { - #include "readTwoPhaseEulerFoamControls.H" - #include "CourantNos.H" + #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 "alphaEqn.H" - #include "liftDragCoeffs.H" - #include "UEqns.H" + #include "UEqn.H" + #include "hEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" - - if (correctAlpha && !pimple.finalIter()) - { - #include "alphaEqn.H" - } } - #include "DDtU.H" - if (pimple.turbCorr()) { - #include "kEpsilon.H" + turbulence->correct(); } } - #include "write.H" + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" 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/Allwclean b/applications/solvers/compressible/rhoSimpleFoam/Allwclean new file mode 100755 index 0000000000..49e4b69482 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/Allwclean @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean +wclean rhoPorousMRFSimpleFoam +wclean rhoSimplecFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/compressible/rhoSimpleFoam/Allwmake b/applications/solvers/compressible/rhoSimpleFoam/Allwmake index 33916c7b39..e0b621c253 100755 --- a/applications/solvers/compressible/rhoSimpleFoam/Allwmake +++ b/applications/solvers/compressible/rhoSimpleFoam/Allwmake @@ -4,5 +4,6 @@ set -x wmake wmake rhoPorousMRFSimpleFoam +wmake rhoSimplecFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H index 21ec2646be..f1bed4d071 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H @@ -3,6 +3,7 @@ tmp UEqn ( fvm::div(phi, U) + - fvm::Sp(fvc::div(phi), U) + turbulence->divDevRhoReff(U) ); 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 8ff402a2b1..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H +++ /dev/null @@ -1,16 +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)") - ); - - hEqn.relax(); - - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 56c444cdab..6c64056be2 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -28,7 +28,7 @@ if (simple.transonic()) ); // Relax the pressure equation to ensure diagonal-dominance - pEqn.relax(mesh.equationRelaxationFactor("pEqn")); + pEqn.relax(); pEqn.setReference(pRefCell, pRefValue); @@ -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/rhoSimplecFoam/Make/files b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/files similarity index 100% rename from applications/solvers/compressible/rhoSimplecFoam/Make/files rename to applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/files diff --git a/applications/solvers/compressible/rhoSimplecFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options similarity index 93% rename from applications/solvers/compressible/rhoSimplecFoam/Make/options rename to applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options index 9d578f011a..30c1e55d31 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/Make/options +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I../rhoSimpleFoam \ + -I.. \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \ diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H similarity index 55% rename from applications/solvers/compressible/rhoSimplecFoam/pEqn.H rename to applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H index a4d9325522..8c7405c346 100644 --- a/applications/solvers/compressible/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,33 +15,33 @@ 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() - + phid*(fvc::interpolate(p) - fvc::interpolate(p, "UD")) - ); - - //refCast(p.boundaryField()[1]).refValue() - // = p.boundaryField()[1]; - fvScalarMatrix pEqn ( fvm::div(phid, p) + fvc::div(phic) - - fvm::Sp(fvc::div(phid), p) - + fvc::div(phid)*p - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); - //pEqn.relax(); + + // Relax the pressure equation to maintain diagonal dominance + pEqn.relax(); pEqn.setReference(pRefCell, pRefValue); @@ -57,22 +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/AU, p) - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); pEqn.setReference(pRefCell, pRefValue); @@ -86,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 @@ -106,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/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C similarity index 95% rename from applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C rename to applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C index 937358e0f5..198b5eea40 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.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 @@ -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/rhoSimplecFoam/createFields.H b/applications/solvers/compressible/rhoSimplecFoam/createFields.H deleted file mode 100644 index 46a382864e..0000000000 --- a/applications/solvers/compressible/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/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/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H index 9db0897efa..06997e52ec 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H @@ -1,18 +1,18 @@ { volScalarField rAU("rAU", 1.0/UEqn.A()); - surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU)); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); + surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, U, phi) - - phig + + phig ); while (pimple.correctNonOrthogonal()) @@ -36,7 +36,7 @@ // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf); + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H index 61d4358a8a..c07a344017 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H @@ -1,12 +1,12 @@ { volScalarField rAU("rAU", 1.0/UEqn().A()); - surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU)); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn().H(); UEqn.clear(); - surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); + surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); surfaceScalarField phiHbyA ( @@ -16,7 +16,7 @@ adjustPhi(phiHbyA, U, p_rgh); - phiHbyA -= phig; + phiHbyA += phig; while (simple.correctNonOrthogonal()) { @@ -39,7 +39,7 @@ // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf); + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index 462b96b17d..96389fee7b 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -6,7 +6,7 @@ thermo.rho() -= psi*p_rgh; volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU)); + surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index d558f6c217..f63e12f363 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -3,13 +3,13 @@ rho.relax(); volScalarField rAU(1.0/UEqn().A()); - surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU)); + surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn().H(); UEqn.clear(); - surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); surfaceScalarField phiHbyA ( @@ -19,7 +19,7 @@ bool closedVolume = adjustPhi(phiHbyA, U, p_rgh); - phiHbyA -= phig; + phiHbyA += phig; while (simple.correctNonOrthogonal()) { @@ -41,7 +41,7 @@ // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf); + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files index 908d7b3a95..522f90e3ef 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files @@ -2,5 +2,4 @@ fluid/compressibleCourantNo.C solid/solidRegionDiffNo.C chtMultiRegionFoam.C -EXE = $(FOAM_USER_APPBIN)/chtMultiRegionFoam - +EXE = $(FOAM_APPBIN)/chtMultiRegionFoam diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H index 73b72af316..d23390edeb 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H @@ -5,7 +5,7 @@ rho.relax(); volScalarField rAU(1.0/UEqn().A()); - surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU)); + surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); U = rAU*UEqn().H(); UEqn.clear(); @@ -15,8 +15,8 @@ dimensionedScalar compressibility = fvc::domainIntegrate(psi); bool compressible = (compressibility.value() > SMALL); - surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - phi -= buoyancyPhi; + surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + phi += phig; // Solve pressure for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) @@ -44,7 +44,7 @@ // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf); + U += rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); U.correctBoundaryConditions(); } } 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/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 76f6f8e30e..f260acc514 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -6,7 +6,7 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn().A()); - surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU)); + surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn().H(); 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/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C index b620536abb..bf12294638 100644 --- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C +++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C @@ -61,6 +61,7 @@ int main(int argc, char *argv[]) fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(fluid.nu(), U) + - (fvc::grad(U) & fvc::grad(fluid.nu())) ); solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H index 2271e85d9b..5a521a04ff 100644 --- a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H @@ -1,5 +1,5 @@ volScalarField rAU(1.0/UEqn().A()); -surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU)); +surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*(UEqn() == sources(U))().H(); diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C index 166afa67b6..557344fa62 100644 --- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C +++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.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 @@ -22,11 +22,11 @@ License along with OpenFOAM. If not, see . Application - uncoupledKinematicParcelFoam + uncoupledKinematicParcelDyMFoam Description Transient solver for the passive transport of a single kinematic - particle could. + particle cloud. Uses a pre-calculated velocity field to evolve the cloud. diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C index 12043167d0..79436625e2 100644 --- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C +++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.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 @@ -26,7 +26,7 @@ Application Description Transient solver for the passive transport of a single kinematic - particle could. + particle cloud. Uses a pre-calculated velocity field to evolve the cloud. diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H index f5364ae314..3f4f485cf6 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H @@ -1,11 +1,11 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); -surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU)); +surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); -surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); +surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); surfaceScalarField phiHbyA ( @@ -15,7 +15,7 @@ surfaceScalarField phiHbyA (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ) - - phig + + phig ); @@ -37,7 +37,7 @@ while (pimple.correctNonOrthogonal()) if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA + p_rghEqn.flux(); - U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/lagrangian/sprayFoam/Allwmake b/applications/solvers/lagrangian/sprayFoam/Allwmake new file mode 100755 index 0000000000..542509db74 --- /dev/null +++ b/applications/solvers/lagrangian/sprayFoam/Allwmake @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake + +wmake sprayEngineFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H index 5f868f3da9..631c9cbba6 100644 --- a/applications/solvers/multiphase/bubbleFoam/pEqn.H +++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H @@ -8,16 +8,21 @@ surfaceScalarField rAU1f(fvc::interpolate(rAU1)); surfaceScalarField rAU2f(fvc::interpolate(rAU2)); - U1 = rAU1*U1Eqn.H(); - U2 = rAU2*U2Eqn.H(); + volVectorField HbyA1("HbyA1", U1); + HbyA1 = rAU1*U1Eqn.H(); + + volVectorField HbyA2("HbyA2", U2); + HbyA2 = rAU2*U2Eqn.H(); surfaceScalarField phiDrag1 ( - fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2 + rAU1f*(g & mesh.Sf()) + fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2 + + rAU1f*(g & mesh.Sf()) ); surfaceScalarField phiDrag2 ( - fvc::interpolate(alpha1/rho2*dragCoef*rAU2)*phi1 + rAU2f*(g & mesh.Sf()) + fvc::interpolate(alpha1/rho2*dragCoef*rAU2)*phi1 + + rAU2f*(g & mesh.Sf()) ); forAll(p.boundaryField(), patchi) @@ -29,16 +34,25 @@ } } - phi1 = (fvc::interpolate(U1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1) - + phiDrag1; - phi2 = (fvc::interpolate(U2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2) - + phiDrag2; + surfaceScalarField phiHbyA1 + ( + (fvc::interpolate(HbyA1) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU1, U1, phi1) + + phiDrag1 + ); - phi = alpha1f*phi1 + alpha2f*phi2; + surfaceScalarField phiHbyA2 + ( + (fvc::interpolate(HbyA2) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU2, U2, phi2) + + phiDrag2 + ); + + surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); surfaceScalarField Dp ( - "(rho*(1|A(U)))", + "Dp", alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2 ); @@ -46,7 +60,7 @@ { fvScalarMatrix pEqn ( - fvm::laplacian(Dp, p) == fvc::div(phi) + fvm::laplacian(Dp, p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); @@ -57,19 +71,17 @@ { surfaceScalarField SfGradp(pEqn.flux()/Dp); - phi1 -= rAU1f*SfGradp/rho1; - phi2 -= rAU2f*SfGradp/rho2; + phi1 = phiHbyA1 - rAU1f*SfGradp/rho1; + phi2 = phiHbyA2 - rAU2f*SfGradp/rho2; phi = alpha1f*phi1 + alpha2f*phi2; p.relax(); SfGradp = pEqn.flux()/Dp; - U1 += (fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1)); - //U1 += rAU1*(fvc::reconstruct(phiDrag1/rAU1f - SfGradp/rho1)); + U1 = HbyA1 + (fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1)); U1.correctBoundaryConditions(); - U2 += (fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2)); - //U2 += rAU2*(fvc::reconstruct(phiDrag2/rAU2f - SfGradp/rho2)); + U2 = HbyA2 + (fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2)); U2.correctBoundaryConditions(); U = alpha1*U1 + alpha2*U2; diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H index 49c7de1473..dbacf1dbd3 100644 --- a/applications/solvers/multiphase/cavitatingFoam/createFields.H +++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H @@ -25,38 +25,6 @@ mesh ); - volScalarField gamma - ( - IOobject - ( - "gamma", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0)) - ); - gamma.oldTime(); - - Info<< "Creating compressibilityModel\n" << endl; - autoPtr psiModel = - barotropicCompressibilityModel::New - ( - thermodynamicProperties, - gamma - ); - - const volScalarField& psi = psiModel->psi(); - - rho == max - ( - psi*p - + (1.0 - gamma)*rhol0 - + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, - rhoMin - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -78,6 +46,27 @@ twoPhaseMixture twoPhaseProperties(U, phiv, "gamma"); + volScalarField& gamma(twoPhaseProperties.alpha1()); + gamma.oldTime(); + + Info<< "Creating compressibilityModel\n" << endl; + autoPtr psiModel = + barotropicCompressibilityModel::New + ( + thermodynamicProperties, + gamma + ); + + const volScalarField& psi = psiModel->psi(); + + rho == max + ( + psi*p + + (1.0 - gamma)*rhol0 + + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, + rhoMin + ); + // Create incompressible turbulence model autoPtr turbulence ( diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H index 34898172ce..ea35f79fa1 100644 --- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H @@ -12,7 +12,7 @@ surfaceScalarField rhof("rhof", fvc::interpolate(rho)); volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU)); + surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); diff --git a/applications/solvers/multiphase/compressibleInterFoam/Allwclean b/applications/solvers/multiphase/compressibleInterFoam/Allwclean index 2b936f934b..2f4544cb4c 100755 --- a/applications/solvers/multiphase/compressibleInterFoam/Allwclean +++ b/applications/solvers/multiphase/compressibleInterFoam/Allwclean @@ -2,6 +2,7 @@ cd ${0%/*} || exit 1 # run from this directory set -x +wclean libso phaseEquationsOfState wclean wclean compressibleInterDyMFoam diff --git a/applications/solvers/multiphase/compressibleInterFoam/Allwmake b/applications/solvers/multiphase/compressibleInterFoam/Allwmake index 644094d070..b4b7f6ffa7 100755 --- a/applications/solvers/multiphase/compressibleInterFoam/Allwmake +++ b/applications/solvers/multiphase/compressibleInterFoam/Allwmake @@ -2,6 +2,7 @@ cd ${0%/*} || exit 1 # run from this directory set -x +wmake libso phaseEquationsOfState wmake wmake compressibleInterDyMFoam diff --git a/applications/solvers/multiphase/compressibleInterFoam/Make/files b/applications/solvers/multiphase/compressibleInterFoam/Make/files index de5437219c..0e009f09bc 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/Make/files +++ b/applications/solvers/multiphase/compressibleInterFoam/Make/files @@ -1,3 +1,4 @@ +derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C compressibleInterFoam.C EXE = $(FOAM_APPBIN)/compressibleInterFoam diff --git a/applications/solvers/multiphase/compressibleInterFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/Make/options index c8ce69c074..ca9a90cf77 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/Make/options +++ b/applications/solvers/multiphase/compressibleInterFoam/Make/options @@ -2,12 +2,14 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -IphaseEquationsOfState/lnInclude \ -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -ltwoPhaseInterfaceProperties \ -lincompressibleTransportModels \ + -lphaseEquationsOfState \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleLESModels \ diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H new file mode 100644 index 0000000000..2605ce345a --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H @@ -0,0 +1,20 @@ +{ + volScalarField kByCv + ( + "kByCv", + (alpha1*k1/Cv1 + alpha2*k2/Cv2) + + (alpha1*rho1 + alpha2*rho2)*turbulence->nut() + ); + + solve + ( + fvm::ddt(rho, T) + + fvm::div(rhoPhi, T) + - fvm::laplacian(kByCv, T) + + p*fvc::div(phi)*(alpha1/Cv1 + alpha2/Cv2) + ); + + // Update compressibilities + psi1 = eos1->psi(p, T); + psi2 = eos2->psi(p, T); +} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options index b9b0ec54da..db259cfbc0 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options @@ -3,6 +3,7 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I../phaseEquationsOfState/lnInclude \ -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ @@ -12,6 +13,7 @@ EXE_INC = \ EXE_LIBS = \ -ltwoPhaseInterfaceProperties \ -lincompressibleTransportModels \ + -lphaseEquationsOfState \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleLESModels \ diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C index 1ac1596c4d..b9f52944e8 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.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 @@ -25,7 +25,7 @@ Application compressibleInterDyMFoam Description - Solver for 2 compressible, isothermal immiscible fluids using a VOF + Solver for 2 compressible, non-isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach, with optional mesh motion and mesh topology changes including adaptive re-meshing. @@ -43,6 +43,7 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" +#include "phaseEquationOfState.H" #include "turbulenceModel.H" #include "pimpleControl.H" @@ -124,11 +125,15 @@ int main(int argc, char *argv[]) solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" + #include "TEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); } } diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H deleted file mode 100644 index 26666c4120..0000000000 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ /dev/null @@ -1,88 +0,0 @@ -{ - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); - - tmp p_rghEqnComp; - - if (pimple.transonic()) - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvm::div(phi, p_rgh) - - fvm::Sp(fvc::div(phi), p_rgh) - ); - } - else - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvc::div(phi, p_rgh) - - fvc::Sp(fvc::div(phi), p_rgh) - ); - } - - - U = rAU*UEqn.H(); - - surfaceScalarField phiU - ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - - phi = phiU + - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix p_rghEqnIncomp - ( - fvc::div(phi) - - fvm::laplacian(rAUf, p_rgh) - ); - - solve - ( - ( - max(alpha1, scalar(0))*(psi1/rho1) - + max(alpha2, scalar(0))*(psi2/rho2) - ) - *p_rghEqnComp() - + p_rghEqnIncomp, - mesh.solver(p_rgh.select(pimple.finalInnerIter())) - ); - - if (pimple.finalNonOrthogonalIter()) - { - dgdt = - (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) - *(p_rghEqnComp & p_rgh); - phi += p_rghEqnIncomp.flux(); - } - } - - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - - p = max - ( - (p_rgh + gh*(alpha1*rho10 + alpha2*rho20)) - /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)), - pMin - ); - - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; - - Info<< "max(U) " << max(mag(U)).value() << endl; - Info<< "min(p_rgh) " << min(p_rgh).value() << endl; - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); -} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index 4ad1b3d01d..7c24e05a1a 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.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 @@ -25,7 +25,7 @@ Application compressibleInterFoam Description - Solver for 2 compressible, isothermal immiscible fluids using a VOF + Solver for 2 compressible, non-isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the "mixture" and a single @@ -40,6 +40,7 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" +#include "phaseEquationOfState.H" #include "turbulenceModel.H" #include "pimpleControl.H" @@ -82,6 +83,7 @@ int main(int argc, char *argv[]) solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" + #include "TEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H index c598cb75ce..1c22600170 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H @@ -12,23 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< "Calculating field alpha1\n" << endl; - volScalarField alpha2("alpha2", scalar(1) - alpha1); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -45,48 +28,20 @@ #include "createPhi.H" - - Info<< "Reading transportProperties\n" << endl; - twoPhaseMixture twoPhaseProperties(U, phi); - - dimensionedScalar rho10 + Info<< "Reading field T\n" << endl; + volScalarField T ( - twoPhaseProperties.subDict + IOobject ( - twoPhaseProperties.phase1Name() - ).lookup("rho0") + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh ); - dimensionedScalar rho20 - ( - twoPhaseProperties.subDict - ( - twoPhaseProperties.phase2Name() - ).lookup("rho0") - ); - - dimensionedScalar psi1 - ( - twoPhaseProperties.subDict - ( - twoPhaseProperties.phase1Name() - ).lookup("psi") - ); - - dimensionedScalar psi2 - ( - twoPhaseProperties.subDict - ( - twoPhaseProperties.phase2Name() - ).lookup("psi") - ); - - dimensionedScalar pMin(twoPhaseProperties.lookup("pMin")); - - Info<< "Calculating field g.h\n" << endl; - volScalarField gh("gh", g & mesh.C()); - surfaceScalarField ghf("ghf", g & mesh.Cf()); - volScalarField p ( IOobject @@ -94,19 +49,120 @@ "p", runTime.timeName(), mesh, - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - max + p_rgh + ); + + + Info<< "Reading transportProperties\n" << endl; + twoPhaseMixture twoPhaseProperties(U, phi); + + volScalarField& alpha1(twoPhaseProperties.alpha1()); + + Info<< "Calculating phase-fraction alpha" << twoPhaseProperties.phase2Name() + << nl << endl; + volScalarField alpha2 + ( + "alpha" + twoPhaseProperties.phase2Name(), + scalar(1) - alpha1 + ); + + dimensionedScalar k1 + ( + "k", + dimensionSet(1, 1, -3, -1, 0), + twoPhaseProperties.subDict ( - (p_rgh + gh*(alpha1*rho10 + alpha2*rho20)) - /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)), - pMin + twoPhaseProperties.phase1Name() + ).lookup("k") + ); + + dimensionedScalar k2 + ( + "k", + dimensionSet(1, 1, -3, -1, 0), + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("k") + ); + + dimensionedScalar Cv1 + ( + "Cv", + dimensionSet(0, 2, -2, -1, 0), + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ).lookup("Cv") + ); + + dimensionedScalar Cv2 + ( + "Cv", + dimensionSet(0, 2, -2, -1, 0), + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("Cv") + ); + + autoPtr eos1 + ( + phaseEquationOfState::New + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ) ) ); - volScalarField rho1(rho10 + psi1*p); - volScalarField rho2(rho20 + psi2*p); + autoPtr eos2 + ( + phaseEquationOfState::New + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ) + ) + ); + + volScalarField psi1 + ( + IOobject + ( + "psi1", + runTime.timeName(), + mesh + ), + eos1->psi(p, T) + ); + psi1.oldTime(); + + volScalarField psi2 + ( + IOobject + ( + "psi2", + runTime.timeName(), + mesh + ), + eos2->psi(p, T) + ); + psi2.oldTime(); + + dimensionedScalar pMin(twoPhaseProperties.lookup("pMin")); + + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + volScalarField rho1("rho1", eos1->rho(p, T)); + volScalarField rho2("rho2", eos2->rho(p, T)); volScalarField rho ( diff --git a/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C new file mode 100644 index 0000000000..e6782e8b3a --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "wallHeatTransferFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(p, iF), + Tinf_(p.size(), 0.0), + alphaWall_(p.size(), 0.0) +{ + refValue() = 0.0; + refGrad() = 0.0; + valueFraction() = 0.0; +} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + mixedFvPatchScalarField(ptf, p, iF, mapper), + Tinf_(ptf.Tinf_, mapper), + alphaWall_(ptf.alphaWall_, mapper) +{} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + mixedFvPatchScalarField(p, iF), + Tinf_("Tinf", dict, p.size()), + alphaWall_("alphaWall", dict, p.size()) +{ + refValue() = Tinf_; + refGrad() = 0.0; + valueFraction() = 0.0; + + if (dict.found("value")) + { + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); + } + else + { + evaluate(); + } +} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& tppsf +) +: + mixedFvPatchScalarField(tppsf), + Tinf_(tppsf.Tinf_), + alphaWall_(tppsf.alphaWall_) +{} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(tppsf, iF), + Tinf_(tppsf.Tinf_), + alphaWall_(tppsf.alphaWall_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::wallHeatTransferFvPatchScalarField::autoMap +( + const fvPatchFieldMapper& m +) +{ + scalarField::autoMap(m); + Tinf_.autoMap(m); + alphaWall_.autoMap(m); +} + + +void Foam::wallHeatTransferFvPatchScalarField::rmap +( + const fvPatchScalarField& ptf, + const labelList& addr +) +{ + mixedFvPatchScalarField::rmap(ptf, addr); + + const wallHeatTransferFvPatchScalarField& tiptf = + refCast(ptf); + + Tinf_.rmap(tiptf.Tinf_, addr); + alphaWall_.rmap(tiptf.alphaWall_, addr); +} + + +void Foam::wallHeatTransferFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const fvPatchScalarField& Cpw = + patch().lookupPatchField("Cp"); + + const fvPatchScalarField& kByCpw = + patch().lookupPatchField("kByCp"); + + valueFraction() = + 1.0/ + ( + 1.0 + + Cpw*kByCpw*patch().deltaCoeffs()/alphaWall_ + ); + + mixedFvPatchScalarField::updateCoeffs(); +} + + +void Foam::wallHeatTransferFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + Tinf_.writeEntry("Tinf", os); + alphaWall_.writeEntry("alphaWall", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField(fvPatchScalarField, wallHeatTransferFvPatchScalarField); +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H new file mode 100644 index 0000000000..cf5284a087 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H @@ -0,0 +1,194 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::wallHeatTransferFvPatchScalarField + +Description + Enthalpy boundary conditions for wall heat transfer + +SourceFiles + wallHeatTransferFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef wallHeatTransferFvPatchScalarField_H +#define wallHeatTransferFvPatchScalarField_H + +#include "mixedFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class wallHeatTransferFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class wallHeatTransferFvPatchScalarField +: + public mixedFvPatchScalarField +{ + // Private data + + //- Tinf + scalarField Tinf_; + + //- alphaWall + scalarField alphaWall_; + + +public: + + //- Runtime type information + TypeName("wallHeatTransfer"); + + + // Constructors + + //- Construct from patch and internal field + wallHeatTransferFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + wallHeatTransferFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given wallHeatTransferFvPatchScalarField + // onto a new patch + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new wallHeatTransferFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new wallHeatTransferFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Access + + //- Return Tinf + const scalarField& Tinf() const + { + return Tinf_; + } + + //- Return reference to Tinf to allow adjustment + scalarField& Tinf() + { + return Tinf_; + } + + //- Return alphaWall + const scalarField& alphaWall() const + { + return alphaWall_; + } + + //- Return reference to alphaWall to allow adjustment + scalarField& alphaWall() + { + return alphaWall_; + } + + + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchScalarField&, + const labelList& + ); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 035e8e237d..9e67a47c4f 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -1,84 +1,97 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + rho1 = eos1->rho(p, T); + rho2 = eos2->rho(p, T); - tmp p_rghEqnComp; + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - if (pimple.transonic()) + tmp p_rghEqnComp1; + tmp p_rghEqnComp2; + + //if (transonic) + //{ + //} + //else { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvm::div(phi, p_rgh) - - fvm::Sp(fvc::div(phi), p_rgh) - ); - } - else - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvc::div(phi, p_rgh) - - fvc::Sp(fvc::div(phi), p_rgh) - ); + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); + surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); + + p_rghEqnComp1 = + fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh)) + + fvc::div(phid1, p_rgh) + - fvc::Sp(fvc::div(phid1), p_rgh); + + p_rghEqnComp2 = + fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh)) + + fvc::div(phid2, p_rgh) + - fvc::Sp(fvc::div(phid2), p_rgh); } + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - U = rAU*UEqn.H(); - - surfaceScalarField phiU + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - phi = phiU + + surfaceScalarField phig + ( ( fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + //thermo.rho() -= psi*p_rgh; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqnIncomp ( - fvc::div(phi) + fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh) ); solve ( ( - max(alpha1, scalar(0))*(psi1/rho1) - + max(alpha2, scalar(0))*(psi2/rho2) + (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1() + + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2() ) - *p_rghEqnComp() + p_rghEqnIncomp, mesh.solver(p_rgh.select(pimple.finalInnerIter())) ); if (pimple.finalNonOrthogonalIter()) { + // Second part of thermodynamic density update + //thermo.rho() += psi*p_rgh; + dgdt = - (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) - *(p_rghEqnComp & p_rgh); - phi += p_rghEqnIncomp.flux(); + ( + pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 + - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 + ); + + phi = phiHbyA + p_rghEqnIncomp.flux(); + + U = HbyA + + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); + p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); - p = max - ( - (p_rgh + gh*(alpha1*rho10 + alpha2*rho20)) - /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)), - pMin - ); - - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; + rho1 = eos1->rho(p, T); + rho2 = eos2->rho(p, T); Info<< "max(U) " << max(mag(U)).value() << endl; Info<< "min(p_rgh) " << min(p_rgh).value() << endl; diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/files b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/files new file mode 100644 index 0000000000..e6e260c74a --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/files @@ -0,0 +1,8 @@ +phaseEquationOfState/phaseEquationOfState.C +phaseEquationOfState/newPhaseEquationOfState.C +constant/constant.C +linear/linear.C +perfectFluid/perfectFluid.C +adiabaticPerfectFluid/adiabaticPerfectFluid.C + +LIB = $(FOAM_LIBBIN)/libphaseEquationsOfState diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/options b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/options new file mode 100644 index 0000000000..0ec1139209 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/Make/options @@ -0,0 +1,6 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude + +LIB_LIBS = \ + -lincompressibleTransportModels diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.C new file mode 100644 index 0000000000..82a195dca9 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.C @@ -0,0 +1,124 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "adiabaticPerfectFluid.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + defineTypeNameAndDebug(adiabaticPerfectFluid, 0); + + addToRunTimeSelectionTable + ( + phaseEquationOfState, + adiabaticPerfectFluid, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::adiabaticPerfectFluid::adiabaticPerfectFluid +( + const dictionary& dict +) +: + phaseEquationOfState(dict), + p0_("p0", dimPressure, dict.lookup("p0")), + rho0_("rho0", dimDensity, dict.lookup("rho0")), + gamma_("gamma", dimless, dict.lookup("gamma")), + B_("B", dimPressure, dict.lookup("B")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::adiabaticPerfectFluid::~adiabaticPerfectFluid() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::phaseEquationsOfState::adiabaticPerfectFluid::rho +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "rho", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + rho0_*pow((p + B_)/(p0_ + B_), 1.0/gamma_) + ) + ); +} + + +Foam::tmp +Foam::phaseEquationsOfState::adiabaticPerfectFluid::psi +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "psi", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + (rho0_/(gamma_*(p0_ + B_))) + *pow((p + B_)/(p0_ + B_), 1.0/gamma_ - 1.0) + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H new file mode 100644 index 0000000000..49f5218e49 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H @@ -0,0 +1,115 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::phaseEquationsOfState::adiabaticPerfectFluid + +Description + AdiabaticPerfectFluid phase density model. + +SourceFiles + adiabaticPerfectFluid.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adiabaticPerfectFluid_H +#define adiabaticPerfectFluid_H + +#include "phaseEquationOfState.H" +#include "dimensionedTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + +/*---------------------------------------------------------------------------*\ + Class adiabaticPerfectFluid Declaration +\*---------------------------------------------------------------------------*/ + +class adiabaticPerfectFluid +: + public phaseEquationOfState +{ + // Private data + + //- Reference pressure + dimensionedScalar p0_; + + //- Reference density + dimensionedScalar rho0_; + + //- The isentropic exponent + dimensionedScalar gamma_; + + //- Pressure offset for a stiffened gas + dimensionedScalar B_; + + +public: + + //- Runtime type information + TypeName("adiabaticPerfectFluid"); + + + // Constructors + + //- Construct from components + adiabaticPerfectFluid + ( + const dictionary& dict + ); + + + //- Destructor + virtual ~adiabaticPerfectFluid(); + + + // Member Functions + + tmp rho + ( + const volScalarField& p, + const volScalarField& T + ) const; + + tmp psi + ( + const volScalarField& p, + const volScalarField& T + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace phaseEquationsOfState +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.C new file mode 100644 index 0000000000..54b6705dd9 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "constant.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + defineTypeNameAndDebug(constant, 0); + + addToRunTimeSelectionTable + ( + phaseEquationOfState, + constant, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::constant::constant +( + const dictionary& dict +) +: + phaseEquationOfState(dict), + rho_("rho", dimDensity, dict.lookup("rho")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::constant::~constant() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::phaseEquationsOfState::constant::rho +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "rho", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + p.mesh(), + rho_ + ) + ); +} + + +Foam::tmp Foam::phaseEquationsOfState::constant::psi +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "psi", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + p.mesh(), + dimensionedScalar("psi", dimDensity/dimPressure, 0) + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.H b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.H new file mode 100644 index 0000000000..5cfe44c26e --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/constant/constant.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::phaseEquationsOfState::constant + +Description + Constant phase density model. + +SourceFiles + constant.C + +\*---------------------------------------------------------------------------*/ + +#ifndef constant_H +#define constant_H + +#include "phaseEquationOfState.H" +#include "dimensionedTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + +/*---------------------------------------------------------------------------*\ + Class constant Declaration +\*---------------------------------------------------------------------------*/ + +class constant +: + public phaseEquationOfState +{ + // Private data + + //- The constant density of the phase + dimensionedScalar rho_; + + +public: + + //- Runtime type information + TypeName("constant"); + + + // Constructors + + //- Construct from components + constant + ( + const dictionary& dict + ); + + + //- Destructor + virtual ~constant(); + + + // Member Functions + + tmp rho + ( + const volScalarField& p, + const volScalarField& T + ) const; + + tmp psi + ( + const volScalarField& p, + const volScalarField& T + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace phaseEquationsOfState +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.C new file mode 100644 index 0000000000..3680735a1b --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "linear.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + defineTypeNameAndDebug(linear, 0); + + addToRunTimeSelectionTable + ( + phaseEquationOfState, + linear, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::linear::linear +( + const dictionary& dict +) +: + phaseEquationOfState(dict), + rho0_("rho0", dimDensity, dict.lookup("rho0")), + psi_("psi", dimDensity/dimPressure, dict.lookup("psi")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::linear::~linear() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::phaseEquationsOfState::linear::rho +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "rho", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + rho0_ + psi_*p + ) + ); +} + + +Foam::tmp Foam::phaseEquationsOfState::linear::psi +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "psi", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + p.mesh(), + psi_ + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.H b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.H new file mode 100644 index 0000000000..d357bb3ec0 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/linear/linear.H @@ -0,0 +1,109 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::phaseEquationsOfState::linear + +Description + Linear phase density model. + +SourceFiles + linear.C + +\*---------------------------------------------------------------------------*/ + +#ifndef linear_H +#define linear_H + +#include "phaseEquationOfState.H" +#include "dimensionedTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + +/*---------------------------------------------------------------------------*\ + Class linear Declaration +\*---------------------------------------------------------------------------*/ + +class linear +: + public phaseEquationOfState +{ + // Private data + + //- The reference density of the phase + dimensionedScalar rho0_; + + //- The constant compressibility of the phase + dimensionedScalar psi_; + + +public: + + //- Runtime type information + TypeName("linear"); + + + // Constructors + + //- Construct from components + linear + ( + const dictionary& dict + ); + + + //- Destructor + virtual ~linear(); + + + // Member Functions + + tmp rho + ( + const volScalarField& p, + const volScalarField& T + ) const; + + tmp psi + ( + const volScalarField& p, + const volScalarField& T + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace phaseEquationsOfState +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.C new file mode 100644 index 0000000000..6a148d8d68 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.C @@ -0,0 +1,119 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "perfectFluid.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + defineTypeNameAndDebug(perfectFluid, 0); + + addToRunTimeSelectionTable + ( + phaseEquationOfState, + perfectFluid, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::perfectFluid::perfectFluid +( + const dictionary& dict +) +: + phaseEquationOfState(dict), + rho0_("rho0", dimDensity, dict.lookup("rho0")), + R_("R", dimensionSet(0, 2, -2, -1, 0), dict.lookup("R")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::phaseEquationsOfState::perfectFluid::~perfectFluid() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::phaseEquationsOfState::perfectFluid::rho +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "rho", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + rho0_ + psi(p, T)*p + ) + ); +} + + +Foam::tmp Foam::phaseEquationsOfState::perfectFluid::psi +( + const volScalarField& p, + const volScalarField& T +) const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "psi", + p.time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + 1.0/(R_*T) + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.H b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.H new file mode 100644 index 0000000000..b854f1d84f --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/perfectFluid/perfectFluid.H @@ -0,0 +1,109 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::phaseEquationsOfState::perfectFluid + +Description + PerfectFluid phase density model. + +SourceFiles + perfectFluid.C + +\*---------------------------------------------------------------------------*/ + +#ifndef perfectFluid_H +#define perfectFluid_H + +#include "phaseEquationOfState.H" +#include "dimensionedTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace phaseEquationsOfState +{ + +/*---------------------------------------------------------------------------*\ + Class perfectFluid Declaration +\*---------------------------------------------------------------------------*/ + +class perfectFluid +: + public phaseEquationOfState +{ + // Private data + + //- The reference density of the phase + dimensionedScalar rho0_; + + //- The fluid constant of the phase + dimensionedScalar R_; + + +public: + + //- Runtime type information + TypeName("perfectFluid"); + + + // Constructors + + //- Construct from components + perfectFluid + ( + const dictionary& dict + ); + + + //- Destructor + virtual ~perfectFluid(); + + + // Member Functions + + tmp rho + ( + const volScalarField& p, + const volScalarField& T + ) const; + + tmp psi + ( + const volScalarField& p, + const volScalarField& T + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace phaseEquationsOfState +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/newPhaseEquationOfState.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/newPhaseEquationOfState.C new file mode 100644 index 0000000000..3d9a842a9f --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/newPhaseEquationOfState.C @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "phaseEquationOfState.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::autoPtr Foam::phaseEquationOfState::New +( + const dictionary& dict +) +{ + word phaseEquationOfStateType + ( + dict.subDict("equationOfState").lookup("type") + ); + + Info<< "Selecting phaseEquationOfState " + << phaseEquationOfStateType << endl; + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(phaseEquationOfStateType); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn("phaseEquationOfState::New") + << "Unknown phaseEquationOfStateType type " + << phaseEquationOfStateType << endl << endl + << "Valid phaseEquationOfState types are : " << endl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return cstrIter()(dict.subDict("equationOfState")); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.C b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.C new file mode 100644 index 0000000000..41ed49322b --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.C @@ -0,0 +1,54 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "phaseEquationOfState.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(phaseEquationOfState, 0); + defineRunTimeSelectionTable(phaseEquationOfState, dictionary); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseEquationOfState::phaseEquationOfState +( + const dictionary& dict +) +: + dict_(dict) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::phaseEquationOfState::~phaseEquationOfState() +{} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.H b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.H new file mode 100644 index 0000000000..45a5079d0d --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/phaseEquationsOfState/phaseEquationOfState/phaseEquationOfState.H @@ -0,0 +1,127 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::phaseEquationOfState + +Description + A2stract base-class for dispersed-phase particle diameter models. + +SourceFiles + phaseEquationOfState.C + newDiameterModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef phaseEquationOfState_H +#define phaseEquationOfState_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "dictionary.H" +#include "volFieldsFwd.H" +#include "runTimeSelectionTables.H" + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class phaseEquationOfState Declaration +\*---------------------------------------------------------------------------*/ + +class phaseEquationOfState +{ +protected: + + // Protected data + + const dictionary& dict_; + + +public: + + //- Runtime type information + TypeName("phaseEquationOfState"); + + + // Declare runtime construction + + declareRunTimeSelectionTable + ( + autoPtr, + phaseEquationOfState, + dictionary, + ( + const dictionary& dict + ), + (dict) + ); + + + // Constructors + + phaseEquationOfState + ( + const dictionary& dict + ); + + + //- Destructor + virtual ~phaseEquationOfState(); + + + // Selectors + + static autoPtr New + ( + const dictionary& dict + ); + + + // Member Functions + + //- Return the phase density + virtual tmp rho + ( + const volScalarField& p, + const volScalarField& T + ) const = 0; + + //- Return the phase compressibility + virtual tmp psi + ( + const volScalarField& p, + const volScalarField& T + ) const = 0; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/DDtU.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/DDtU.H index 96b6a63a43..1685caa588 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/DDtU.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/DDtU.H @@ -3,9 +3,11 @@ fvc::ddt(U1) + fvc::div(phi1, U1) - fvc::div(phi1)*U1; + mrfZones.addCoriolis(U1, DDtU1); DDtU2 = fvc::ddt(U2) + fvc::div(phi2, U2) - fvc::div(phi2)*U2; + mrfZones.addCoriolis(U2, DDtU2); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H index ed150abb88..5d61f0eaf7 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H @@ -18,7 +18,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); volTensorField Rc1 ( "Rc1", - ((2.0/3.0)*I)*(sqr(Ct)*k + nuEff1*tr(gradU1T)) - nuEff1*gradU1T + (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T ); if (kineticTheory.on()) @@ -28,19 +28,24 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); U1Eqn = ( - (scalar(1) + Cvm*rho2*alpha2/rho1)* + fvm::ddt(alpha1, U1) + + fvm::div(alphaPhi1, U1) + - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1) + + + Cvm*rho2*alpha1*alpha2/rho1* ( - fvm::ddt(alpha1, U1) - + fvm::div(alphaPhi1, U1) - - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1) + fvm::ddt(U1) + + fvm::div(phi1, U1) + - fvm::Sp(fvc::div(phi1), U1) ) + - fvm::laplacian(alpha1*nuEff1, U1) + fvc::div(alpha1*Rc1) == - fvm::Sp(dragCoeff/rho1, U1) - alpha1*alpha2/rho1*(liftForce - Cvm*rho2*DDtU2) ); - + mrfZones.addCoriolis(alpha1*(1 + Cvm*rho2*alpha2/rho1), U1Eqn); U1Eqn.relax(); } @@ -49,24 +54,29 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); volTensorField Rc2 ( "Rc2", - ((2.0/3.0)*I)*(k + nuEff2*tr(gradU2T)) - nuEff2*gradU2T + (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T ); U2Eqn = ( - (scalar(1) + Cvm*rho2*alpha1/rho2)* + fvm::ddt(alpha2, U2) + + fvm::div(alphaPhi2, U2) + - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2) + + + Cvm*rho2*alpha1*alpha2/rho2* ( - fvm::ddt(alpha2, U2) - + fvm::div(alphaPhi2, U2) - - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2) + fvm::ddt(U2) + + fvm::div(phi2, U2) + - fvm::Sp(fvc::div(phi2), U2) ) + - fvm::laplacian(alpha2*nuEff2, U2) + fvc::div(alpha2*Rc2) == - fvm::Sp(dragCoeff/rho2, U2) + alpha1*alpha2/rho2*(liftForce + Cvm*rho2*DDtU1) ); - + mrfZones.addCoriolis(alpha2*(1 + Cvm*rho2*alpha1/rho2), U2Eqn); U2Eqn.relax(); } } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C index d29da482ba..c736f8194d 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.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 @@ -33,6 +33,7 @@ Description #include "fvCFD.H" #include "nearWallDist.H" #include "wallFvPatch.H" +#include "fixedValueFvsPatchFields.H" #include "Switch.H" #include "IFstream.H" @@ -45,6 +46,8 @@ Description #include "pimpleControl.H" +#include "MRFZones.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) @@ -55,6 +58,7 @@ int main(int argc, char *argv[]) #include "createMesh.H" #include "readGravitationalAcceleration.H" #include "createFields.H" + #include "createMRFZones.H" #include "readPPProperties.H" #include "initContinuityErrs.H" #include "readTimeControls.H" diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H index 93c35ee1e1..c773650be2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H @@ -310,14 +310,21 @@ Info << "dispersedPhase is " << dispersedPhase << endl; - scalar minInterfaceAlpha + scalar residualPhaseFraction ( readScalar ( - interfacialProperties.lookup("minInterfaceAlpha") + interfacialProperties.lookup("residualPhaseFraction") ) ); + dimensionedScalar residualSlip + ( + "residualSlip", + dimVelocity, + interfacialProperties.lookup("residualSlip") + ); + kineticTheoryModel kineticTheory ( phase1, diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createMRFZones.H similarity index 73% rename from applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H rename to applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createMRFZones.H index 2fc0eff9eb..4d5c2bab72 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createMRFZones.H @@ -1,3 +1,4 @@ MRFZones mrfZones(mesh); mrfZones.correctBoundaryVelocity(U1); mrfZones.correctBoundaryVelocity(U2); + mrfZones.correctBoundaryVelocity(U); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H index 0f28f2330d..d6ccf90289 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H @@ -36,7 +36,7 @@ volScalarField heatTransferCoeff { volVectorField Ur(U1 - U2); - volScalarField magUr(mag(Ur)); + volScalarField magUr(mag(Ur) + residualSlip); if (dispersedPhase == "1") { @@ -69,12 +69,20 @@ volScalarField heatTransferCoeff << exit(FatalError); } - volScalarField alpha1Coeff - ( - (alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha) - ); - dragCoeff *= alpha1Coeff; - heatTransferCoeff *= alpha1Coeff; + volScalarField alphaCoeff(max(alpha1*alpha2, residualPhaseFraction)); + dragCoeff *= alphaCoeff; + heatTransferCoeff *= alphaCoeff; liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U)); + + // Remove lift, drag and phase heat-transfer at fixed-flux boundaries + forAll(phi1.boundaryField(), patchi) + { + if (isA(phi1.boundaryField()[patchi])) + { + dragCoeff.boundaryField()[patchi] = 0.0; + heatTransferCoeff.boundaryField()[patchi] = 0.0; + liftForce.boundaryField()[patchi] = vector::zero; + } + } } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index d2ef15776d..11f3d4ddb5 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -2,6 +2,11 @@ rho1 = rho10 + psi1*p; rho2 = rho20 + psi2*p; + mrfZones.absoluteFlux(phi1.oldTime()); + mrfZones.absoluteFlux(phi1); + mrfZones.absoluteFlux(phi2.oldTime()); + mrfZones.absoluteFlux(phi2); + surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField alpha2f(scalar(1) - alpha1f); @@ -25,6 +30,7 @@ + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2 + rAlphaAU1f*(g & mesh.Sf()) ); + mrfZones.relativeFlux(phiHbyA1); surfaceScalarField phiHbyA2 ( @@ -34,9 +40,18 @@ + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1 + rAlphaAU2f*(g & mesh.Sf()) ); + mrfZones.relativeFlux(phiHbyA2); + + mrfZones.relativeFlux(phi1.oldTime()); + mrfZones.relativeFlux(phi1); + mrfZones.relativeFlux(phi2.oldTime()); + mrfZones.relativeFlux(phi2); surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); + HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2; + HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1; + surfaceScalarField Dp ( "Dp", @@ -91,14 +106,8 @@ { surfaceScalarField mSfGradp = pEqnIncomp.flux()/Dp; - phi1.boundaryField() == - (fvc::interpolate(U1) & mesh.Sf())().boundaryField(); phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1); - - phi2.boundaryField() == - (fvc::interpolate(U2) & mesh.Sf())().boundaryField(); phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2); - phi = alpha1f*phi1 + alpha2f*phi2; dgdt = @@ -110,23 +119,25 @@ p.relax(); mSfGradp = pEqnIncomp.flux()/Dp; - volVectorField U10 = U1; - U1 = HbyA1 - + (1.0/rho1)*rAU1*dragCoeff*U2 + fvc::reconstruct ( - rAlphaAU1f*(g & mesh.Sf()) - + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1) + rAlphaAU1f + *( + (g & mesh.Sf()) + + mSfGradp/fvc::interpolate(rho1) + ) ); U1.correctBoundaryConditions(); U2 = HbyA2 - + (1.0/rho2)*rAU2*dragCoeff*U10 + fvc::reconstruct ( - rAlphaAU2f*(g & mesh.Sf()) - + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2) + rAlphaAU2f + *( + (g & mesh.Sf()) + + mSfGradp/fvc::interpolate(rho2) + ) ); U2.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H index 6fc51c7cee..9c56c4d26f 100644 --- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H @@ -1,29 +1,35 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); - surfaceScalarField phiU + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - mrfZones.relativeFlux(phiU); + mrfZones.relativeFlux(phiHbyA); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = phiU + + surfaceScalarField phig ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -32,13 +38,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H index 8dbb489c8d..fa56de00bb 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H @@ -1,5 +1,4 @@ label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); - label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); if (nAlphaSubCycles > 1) diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H index 34dea334e8..68765262db 100644 --- a/applications/solvers/multiphase/interFoam/createFields.H +++ b/applications/solvers/multiphase/interFoam/createFields.H @@ -12,20 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -46,6 +32,8 @@ Info<< "Reading transportProperties\n" << endl; twoPhaseMixture twoPhaseProperties(U, phi); + volScalarField& alpha1(twoPhaseProperties.alpha1()); + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H index d94e075424..d942f2c501 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H @@ -1,31 +1,39 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - phiAbs = - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phiAbs); + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phiAbs) + ); if (p_rgh.needReference()) { - fvc::makeRelative(phiAbs, U); - adjustPhi(phiAbs, U, p_rgh); - fvc::makeAbsolute(phiAbs, U); + fvc::makeRelative(phiHbyA, U); + adjustPhi(phiHbyA, U, p_rgh); + fvc::makeAbsolute(phiHbyA, U); } - phi = phiAbs + + surfaceScalarField phig ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -34,13 +42,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiAbs)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" phiAbs = phi; diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H index a88b5627e2..d1bc7bf1ed 100644 --- a/applications/solvers/multiphase/interFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/pEqn.H @@ -1,28 +1,34 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); - surfaceScalarField phiU + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = phiU + + surfaceScalarField phig ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -31,13 +37,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H index c9b65eb9d5..e3ef0dc9ef 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H @@ -1,24 +1,28 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - surfaceScalarField phiU + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = - phiU - + ( + surfaceScalarField phig + ( + ( fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; Pair > vDotP = twoPhaseProperties->vDotP(); const volScalarField& vDotcP = vDotP[0](); @@ -28,7 +32,7 @@ { fvScalarMatrix p_rghEqn ( - fvc::div(phi) - fvm::laplacian(rAUf, p_rgh) + fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh) - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh) ); @@ -38,13 +42,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi += p_rghEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H b/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H index f0f5bad356..dc39b70863 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H @@ -6,4 +6,6 @@ forAllIter(PtrDictionary, fluid.phases(), iter) fvc::ddt(phase.U()) + fvc::div(phase.phi(), phase.U()) - fvc::div(phase.phi())*phase.U(); + + mrfZones.addCoriolis(phase.U(), phase.DDtU()); } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H index bce95f7d5a..d4b371915c 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H @@ -15,12 +15,17 @@ forAllIter(PtrDictionary, fluid.phases(), iter) phasei, new fvVectorMatrix ( - (scalar(1) + fluid.Cvm(phase)/phase.rho())* + fvm::ddt(alpha, U) + + fvm::div(phase.phiAlpha(), U) + - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U) + + + (alpha/phase.rho())*fluid.Cvm(phase)* ( - fvm::ddt(alpha, U) - + fvm::div(phase.phiAlpha(), U) - - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U) + fvm::ddt(U) + + fvm::div(phase.phi(), U) + - fvm::Sp(fvc::div(phase.phi()), U) ) + - fvm::laplacian(alpha*nuEff, U) - fvc::div ( @@ -44,7 +49,11 @@ forAllIter(PtrDictionary, fluid.phases(), iter) ) ) ); - mrfZones.addCoriolis(alpha, UEqns[phasei]); + mrfZones.addCoriolis + ( + alpha*(1 + (1/phase.rho())*fluid.Cvm(phase)), + UEqns[phasei] + ); UEqns[phasei].relax(); phasei++; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H index 78f93384b0..57a79b6091 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H @@ -32,7 +32,7 @@ dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0); - adjustPhi(phi, U, pcorr); + //adjustPhi(phi, U, pcorr); while (pimple.correctNonOrthogonal()) { diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 2224030613..620499a2e7 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.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 @@ -25,6 +25,7 @@ License #include "multiphaseSystem.H" #include "alphaContactAngleFvPatchScalarField.H" +#include "fixedValueFvsPatchFields.H" #include "Time.H" #include "subCycle.H" #include "MULES.H" @@ -610,6 +611,21 @@ Foam::tmp Foam::multiphaseSystem::Svm } } + // Remove lift at fixed-flux boundaries + forAll(phase.phi().boundaryField(), patchi) + { + if + ( + isA + ( + phase.phi().boundaryField()[patchi] + ) + ) + { + tSvm().boundaryField()[patchi] = vector::zero; + } + } + return tSvm; } @@ -623,9 +639,7 @@ Foam::multiphaseSystem::dragCoeffs() const { const dragModel& dm = *iter(); - dragCoeffsPtr().insert - ( - iter.key(), + volScalarField* Kptr = ( max ( @@ -642,8 +656,24 @@ Foam::multiphaseSystem::dragCoeffs() const dm.residualSlip() ) ) - ).ptr() - ); + ).ptr(); + + // Remove drag at fixed-flux boundaries + forAll(dm.phase1().phi().boundaryField(), patchi) + { + if + ( + isA + ( + dm.phase1().phi().boundaryField()[patchi] + ) + ) + { + Kptr->boundaryField()[patchi] = 0.0; + } + } + + dragCoeffsPtr().insert(iter.key(), Kptr); } return dragCoeffsPtr; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index 4a817a1001..6dc67ce353 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -76,6 +76,7 @@ + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi()) + rAlphaAUfs[phasei]*(g & mesh.Sf()) ); + mrfZones.relativeFlux(phiHbyAs[phasei]); multiphaseSystem::dragModelTable::const_iterator dmIter = @@ -89,11 +90,20 @@ ++dmIter, ++dcIter ) { - const phaseModel *phase2Ptr = &dmIter()->phase2(); - if (phase2Ptr == &phase) + const phaseModel *phase2Ptr = NULL; + + if (&phase == &dmIter()->phase1()) + { + phase2Ptr = &dmIter()->phase2(); + } + else if (&phase == &dmIter()->phase2()) { phase2Ptr = &dmIter()->phase1(); } + else + { + continue; + } phiHbyAs[phasei] += fvc::interpolate @@ -105,6 +115,9 @@ (1.0/phase.rho())*rAUs[phasei]*(*dcIter()) *phase2Ptr->U(); + // Alternative flux-pressure consistent drag + // but not momentum conservative + // // HbyAs[phasei] += fvc::reconstruct // ( // fvc::interpolate @@ -119,6 +132,17 @@ phasei++; } + phasei = 0; + forAllIter(PtrDictionary, fluid.phases(), iter) + { + phaseModel& phase = iter(); + + mrfZones.relativeFlux(phase.phi().oldTime()); + mrfZones.relativeFlux(phase.phi()); + + phasei++; + } + surfaceScalarField Dp ( IOobject @@ -135,11 +159,10 @@ forAllIter(PtrDictionary, fluid.phases(), iter) { phaseModel& phase = iter(); - Dp += alphafs[phasei]*rAlphaAUfs[phasei]/phase.rho(); + Dp += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho(); phasei++; } - Dp = mag(Dp); while (pimple.correctNonOrthogonal()) { @@ -149,6 +172,8 @@ - fvm::laplacian(Dp, p) ); + pEqnIncomp.setReference(pRefCell, pRefValue); + solve ( // ( @@ -172,8 +197,10 @@ phase.phi() = phiHbyAs[phasei] + rAlphaAUfs[phasei]*mSfGradp/phase.rho(); - mrfZones.relativeFlux(phase.phi().oldTime()); - phi += alphafs[phasei]*phase.phi(); + phi += + alphafs[phasei]*phiHbyAs[phasei] + + mag(alphafs[phasei]*rAlphaAUfs[phasei]) + *mSfGradp/phase.rho(); phasei++; } diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H index efb181f3f2..7bc1eeb12c 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H @@ -1,31 +1,32 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - surfaceScalarField phiU + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) - //+ fvc::ddtPhiCorr(rAU, rho, U, phi) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - mrfZones.relativeFlux(phiU); + mrfZones.relativeFlux(phiHbyA); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = - phiU - + ( - mixture.surfaceTensionForce() - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + surfaceScalarField phig + ( + - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -34,13 +35,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H index e7e5bc14cf..c1e346ce43 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H @@ -1,29 +1,34 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); - surfaceScalarField phiU + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = phiU + + surfaceScalarField phig ( - mixture.surfaceTensionForce() - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -32,13 +37,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H index d89fd3ef4e..9f27e992a8 100644 --- a/applications/solvers/multiphase/settlingFoam/pEqn.H +++ b/applications/solvers/multiphase/settlingFoam/pEqn.H @@ -1,54 +1,61 @@ -volScalarField rAU(1.0/UEqn.A()); - -surfaceScalarField rAUf -( - "(rho*(1|A(U)))", - fvc::interpolate(rho)*fvc::interpolate(rAU) -); - -U = rAU*UEqn.H(); -phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - -surfaceScalarField phiU("phiU", phi); -phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf(); - -while (pimple.correctNonOrthogonal()) { - fvScalarMatrix p_rghEqn + volScalarField rAU("rAU", 1.0/UEqn.A()); + + surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA ( - fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi) + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) ); - p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + surfaceScalarField phig + ( + - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf() + ); - p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); + phiHbyA += phig; - if (pimple.finalNonOrthogonalIter()) + while (pimple.correctNonOrthogonal()) { - phi -= p_rghEqn.flux(); + fvScalarMatrix p_rghEqn + ( + fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA) + ); + + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); + } } + + p == p_rgh + rho*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; + } + + #include "rhoEqn.H" + #include "compressibleContinuityErrs.H" } - -p == p_rgh + rho*gh; - -if (p_rgh.needReference()) -{ - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, pRefCell) - ); - p_rgh = p - rho*gh; -} - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -U += rAU*fvc::reconstruct((phi - phiU)/rAUf); -U.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H index 73af502ce2..0d01b9a9e5 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H @@ -12,20 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -45,6 +31,8 @@ Info<< "Reading transportProperties\n" << endl; twoPhaseMixture twoPhaseProperties(U, phi); + volScalarField& alpha1(twoPhaseProperties.alpha1()); + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H index ac7fc35f68..ddfca4e3ea 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H @@ -1,24 +1,31 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); - U = rAU*UEqn.H(); - surfaceScalarField phiU + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiU, U, p_rgh); + adjustPhi(phiHbyA, U, p_rgh); - phi = phiU - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf(); + surfaceScalarField phig + ( + - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf() + ); + + phiHbyA += phig; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); @@ -27,13 +34,13 @@ if (pimple.finalNonOrthogonalIter()) { - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - #include "continuityErrs.H" p == p_rgh + rho*gh; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean index 63db39ff05..cc138bc068 100755 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean @@ -6,6 +6,5 @@ wclean libso phaseModel wclean libso interfacialModels wclean libso kineticTheoryModels wclean -wclean MRFtwoPhaseEulerFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake index faf438d0bd..29294d166a 100755 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake @@ -6,6 +6,5 @@ wmake libso phaseModel wmake libso interfacialModels wmake libso kineticTheoryModels wmake -wmake MRFtwoPhaseEulerFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H index 96b6a63a43..1685caa588 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H @@ -3,9 +3,11 @@ fvc::ddt(U1) + fvc::div(phi1, U1) - fvc::div(phi1)*U1; + mrfZones.addCoriolis(U1, DDtU1); DDtU2 = fvc::ddt(U2) + fvc::div(phi2, U2) - fvc::div(phi2)*U2; + mrfZones.addCoriolis(U2, DDtU2); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files deleted file mode 100644 index 45960722ae..0000000000 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -MRFTwoPhaseEulerFoam.C - -EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options deleted file mode 100644 index b9b19059da..0000000000 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options +++ /dev/null @@ -1,17 +0,0 @@ -EXE_INC = \ - -I.. \ - -I../../bubbleFoam \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ - -IturbulenceModel \ - -I../kineticTheoryModels/lnInclude \ - -I../interfacialModels/lnInclude \ - -I../phaseModel/lnInclude \ - -EXE_LIBS = \ - -lEulerianInterfacialModels \ - -lfiniteVolume \ - -lmeshTools \ - -lincompressibleTransportModels \ - -lphaseModel \ - -lkineticTheoryModel diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H index 8d09ccd3ef..6fae832a21 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H @@ -53,7 +53,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); //+ alpha2/rho1*K*U2 // Explicit drag transfered to p-equation - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2) ); - + mrfZones.addCoriolis(scalar(1) + Cvm*rho2*alpha2/rho1, U1Eqn); U1Eqn.relax(); } @@ -93,7 +93,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); //+ alpha1/rho2*K*U1 // Explicit drag transfered to p-equation + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1) ); - + mrfZones.addCoriolis(scalar(1) + Cvm*rho2*alpha1/rho2, U2Eqn); U2Eqn.relax(); } } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H index 52f09e4793..34fdc24afc 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H @@ -1,6 +1,9 @@ { - word scheme("div(phi,alpha1)"); - word schemer("div(phir,alpha1)"); + label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); + label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); + + word alphaScheme("div(phi,alpha1)"); + word alpharScheme("div(phir,alpha1)"); surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); @@ -15,14 +18,39 @@ for (int acorr=0; acorr alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + surfaceScalarField phiAlpha + ( + fvc::flux + ( + phic, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, alpha2, alpharScheme), + alpha1, + alpharScheme + ) + ); - if (g0.value() > 0.0) + MULES::explicitSolve + ( + alpha1, + phi, + phiAlpha, + (g0.value() > 0 ? alphaMax : 1), + 0 + ); + } + + if (g0.value() > 0) { ppMagf = rAU1f*fvc::interpolate ( @@ -30,19 +58,23 @@ *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) ); - alpha1Eqn -= fvm::laplacian + fvScalarMatrix alpha1Eqn ( - (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, - alpha1, - "laplacian(alpha1PpMag,alpha1)" + fvm::ddt(alpha1) - fvc::ddt(alpha1) + - fvm::laplacian + ( + (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, + alpha1, + "laplacian(alpha1PpMag,alpha1)" + ) ); + + alpha1Eqn.relax(); + alpha1Eqn.solve(); + + #include "packingLimiter.H" } - alpha1Eqn.relax(); - alpha1Eqn.solve(); - - #include "packingLimiter.H" - alpha2 = scalar(1) - alpha1; Info<< "Dispersed phase volume fraction = " diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index b71c7849e5..314d9b55ea 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H @@ -93,16 +93,22 @@ dimensionedScalar Cvm ( + "Cvm", + dimless, transportProperties.lookup("Cvm") ); dimensionedScalar Cl ( + "Cl", + dimless, transportProperties.lookup("Cl") ); dimensionedScalar Ct ( + "Ct", + dimless, transportProperties.lookup("Ct") ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H new file mode 100644 index 0000000000..4d5c2bab72 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createMRFZones.H @@ -0,0 +1,4 @@ + MRFZones mrfZones(mesh); + mrfZones.correctBoundaryVelocity(U1); + mrfZones.correctBoundaryVelocity(U2); + mrfZones.correctBoundaryVelocity(U); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H b/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H index 097123378c..415afe9483 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H @@ -19,3 +19,13 @@ ( Cl*(alpha2*rho2 + alpha1*rho1)*(Ur ^ fvc::curl(U)) ); + + // Remove lift and drag at fixed-flux boundaries + forAll(phi1.boundaryField(), patchi) + { + if (isA(phi1.boundaryField()[patchi])) + { + K.boundaryField()[patchi] = 0.0; + liftCoeff.boundaryField()[patchi] = vector::zero; + } + } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H new file mode 100644 index 0000000000..9ce4537979 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H @@ -0,0 +1,11 @@ +{ + DDtU1 = + fvc::ddt(U1) + + fvc::div(phi, U1) + - fvc::div(phi)*U1; + + DDtU2 = + fvc::ddt(U2) + + fvc::div(phi, U2) + - fvc::div(phi)*U2; +} diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H similarity index 72% rename from applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H index 0c0cc1543a..cae326e44a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H @@ -2,6 +2,16 @@ fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime); fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); { + volScalarField Cvma + ( + Cvm + /( + 1 + + mag(fvc::grad(alpha1)) + *dimensionedScalar("l", dimLength, 1e-2) + ) + ); + { volTensorField gradU1T(T(fvc::grad(U1))); @@ -34,11 +44,15 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); U1Eqn = ( - (scalar(1) + Cvm*rho2*alpha2/rho1)* + fvm::ddt(U1) + + fvm::div(phi1, U1) + - fvm::Sp(fvc::div(phi1), U1) + + + Cvma*alpha2*rho/rho1* ( fvm::ddt(U1) - + fvm::div(phi1, U1, "div(phi1,U1)") - - fvm::Sp(fvc::div(phi1), U1) + + fvm::div(phi, U1) + - fvm::Sp(fvc::div(phi), U1) ) - fvm::laplacian(nuEff1, U1) @@ -51,9 +65,11 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); // g // Buoyancy term transfered to p-equation - fvm::Sp(alpha2/rho1*K, U1) //+ alpha2/rho1*K*U2 // Explicit drag transfered to p-equation - - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2) + // - fvm::Sp(K/rho1, U1) + ////+ K/rho1*U // Explicit drag transfered to p-equation + - alpha2/rho1*liftCoeff + Cvma*alpha2*rho*DDtU2/rho1 ); - mrfZones.addCoriolis(U1Eqn); + U1Eqn.relax(); } @@ -73,11 +89,15 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); U2Eqn = ( - (scalar(1) + Cvm*rho2*alpha1/rho2)* + fvm::ddt(U2) + + fvm::div(phi2, U2) + - fvm::Sp(fvc::div(phi2), U2) + + + Cvma*alpha1*rho/rho2* ( fvm::ddt(U2) - + fvm::div(phi2, U2, "div(phi2,U2)") - - fvm::Sp(fvc::div(phi2), U2) + + fvm::div(phi, U2) + - fvm::Sp(fvc::div(phi), U2) ) - fvm::laplacian(nuEff2, U2) @@ -85,15 +105,16 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); + fvm::div(phiR2, U2, "div(phi2,U2)") - fvm::Sp(fvc::div(phiR2), U2) - + (fvc::grad(alpha2)/(fvc::average(alpha2) + scalar(0.001)) & Rc2) == // g // Buoyancy term transfered to p-equation - fvm::Sp(alpha1/rho2*K, U2) //+ alpha1/rho2*K*U1 // Explicit drag transfered to p-equation - + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1) + // - fvm::Sp(K/rho2, U2) + ////+ K/rho2*U // Explicit drag transfered to p-equation + + alpha1/rho2*liftCoeff + Cvma*alpha1*rho*DDtU1/rho2 ); - mrfZones.addCoriolis(U2Eqn); + U2Eqn.relax(); } } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index 91f8302d7f..9b543bcd25 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -5,8 +5,8 @@ volScalarField rAU1(1.0/U1Eqn.A()); volScalarField rAU2(1.0/U2Eqn.A()); - rAU1f = fvc::interpolate(rAU1); - surfaceScalarField rAU2f(fvc::interpolate(rAU2)); + rAU1f = 1.0/fvc::interpolate(U1Eqn.A()); + surfaceScalarField rAU2f(1.0/fvc::interpolate(U2Eqn.A())); volVectorField HbyA1("HbyA1", U1); HbyA1 = rAU1*U1Eqn.H(); @@ -14,34 +14,21 @@ volVectorField HbyA2("HbyA2", U2); HbyA2 = rAU2*U2Eqn.H(); - surfaceScalarField phiDrag1 - ( - fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf()) - ); + mrfZones.absoluteFlux(phi1.oldTime()); + mrfZones.absoluteFlux(phi1); + mrfZones.absoluteFlux(phi2.oldTime()); + mrfZones.absoluteFlux(phi2); + + surfaceScalarField ppDrag("ppDrag", 0.0*phi1); if (g0.value() > 0.0) { - phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf(); + ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf(); } if (kineticTheory.on()) { - phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf(); - } - - surfaceScalarField phiDrag2 - ( - fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf()) - ); - - // Fix for gravity on outlet boundary. - forAll(p.boundaryField(), patchi) - { - if (isA(p.boundaryField()[patchi])) - { - phiDrag1.boundaryField()[patchi] = 0.0; - phiDrag2.boundaryField()[patchi] = 0.0; - } + ppDrag -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf(); } surfaceScalarField phiHbyA1 @@ -49,19 +36,32 @@ "phiHbyA1", (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1) - + phiDrag1 + + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + + ppDrag + + rAU1f*(g & mesh.Sf()) ); + mrfZones.relativeFlux(phiHbyA1); surfaceScalarField phiHbyA2 ( "phiHbyA2", (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2) - + phiDrag2 + + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + + rAU2f*(g & mesh.Sf()) ); + mrfZones.relativeFlux(phiHbyA2); + + mrfZones.relativeFlux(phi1.oldTime()); + mrfZones.relativeFlux(phi1); + mrfZones.relativeFlux(phi2.oldTime()); + mrfZones.relativeFlux(phi2); surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); + HbyA1 += alpha2*(1.0/rho1)*rAU1*K*U2; + HbyA2 += alpha1*(1.0/rho2)*rAU2*K*U1; + surfaceScalarField Dp ( "Dp", @@ -83,23 +83,26 @@ { surfaceScalarField SfGradp(pEqn.flux()/Dp); - phi1.boundaryField() == - (fvc::interpolate(U1) & mesh.Sf())().boundaryField(); phi1 = phiHbyA1 - rAU1f*SfGradp/rho1; - - phi2.boundaryField() == - (fvc::interpolate(U2) & mesh.Sf())().boundaryField(); phi2 = phiHbyA2 - rAU2f*SfGradp/rho2; - phi = alpha1f*phi1 + alpha2f*phi2; p.relax(); SfGradp = pEqn.flux()/Dp; - U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1); + U1 = HbyA1 + + fvc::reconstruct + ( + ppDrag + + rAU1f*((g & mesh.Sf()) - SfGradp/rho1) + ); U1.correctBoundaryConditions(); - U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2); + U2 = HbyA2 + + fvc::reconstruct + ( + rAU2f*((g & mesh.Sf()) - SfGradp/rho2) + ); U2.correctBoundaryConditions(); U = alpha1*U1 + alpha2*U2; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old similarity index 87% rename from applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old index 079ecfefe4..348cc847d7 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old @@ -69,6 +69,11 @@ ); mrfZones.relativeFlux(phiHbyA2); + mrfZones.relativeFlux(phi1.oldTime()); + mrfZones.relativeFlux(phi1); + mrfZones.relativeFlux(phi2.oldTime()); + mrfZones.relativeFlux(phi2); + surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); surfaceScalarField Dp @@ -92,18 +97,8 @@ { surfaceScalarField SfGradp(pEqn.flux()/Dp); - phi1.boundaryField() == - (fvc::interpolate(U1) & mesh.Sf())().boundaryField(); - mrfZones.relativeFlux(phi1); phi1 = phiHbyA1 - rAU1f*SfGradp/rho1; - mrfZones.relativeFlux(phi1.oldTime()); - - phi2.boundaryField() == - (fvc::interpolate(U2) & mesh.Sf())().boundaryField(); - mrfZones.relativeFlux(phi2); phi2 = phiHbyA2 - rAU2f*SfGradp/rho2; - mrfZones.relativeFlux(phi2.oldTime()); - phi = alpha1f*phi1 + alpha2f*phi2; p.relax(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C index d81ad8e703..aebc7ba6cb 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C @@ -43,14 +43,20 @@ Foam::phaseModel::phaseModel name_(phaseName), d_ ( + "d", + dimLength, dict_.lookup("d") ), nu_ ( + "nu", + dimensionSet(0, 2, -1, 0, 0), dict_.lookup("nu") ), rho_ ( + "rho", + dimDensity, dict_.lookup("rho") ), U_ diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H index a345c4e53b..4e6847debe 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H @@ -1,5 +1,3 @@ #include "readTimeControls.H" - int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr"))); - Switch correctAlpha(pimple.dict().lookup("correctAlpha")); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index 19810a15f1..178032db6f 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.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 @@ -31,8 +31,11 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "MULES.H" +#include "subCycle.H" #include "nearWallDist.H" #include "wallFvPatch.H" +#include "fixedValueFvsPatchFields.H" #include "Switch.H" #include "IFstream.H" @@ -43,6 +46,7 @@ Description #include "kineticTheoryModel.H" #include "pimpleControl.H" +#include "MRFZones.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -56,6 +60,7 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "readPPProperties.H" #include "initContinuityErrs.H" + #include "createMRFZones.H" #include "readTimeControls.H" #include "CourantNo.H" #include "setInitialDeltaT.H" diff --git a/applications/test/LduMatrix/LduMatrixTest.C b/applications/test/LduMatrix/LduMatrixTest.C new file mode 100644 index 0000000000..13038bef45 --- /dev/null +++ b/applications/test/LduMatrix/LduMatrixTest.C @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 . + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "fvMesh.H" +#include "volFields.H" +#include "LduMatrix.H" +#include "diagTensorField.H" +#include "TPCG.H" +#include "TPBiCG.H" +#include "NoPreconditioner.H" + +using namespace Foam; + +typedef Foam::LduMatrix + lduVectorMatrix; +defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0); + +typedef Foam::DiagonalSolver + lduVectorDiagonalSolver; +defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0); + +template<> +const vector lduVectorMatrix::great_(1e15, 1e15, 1e15); + +template<> +const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15); + +namespace Foam +{ + typedef LduMatrix::preconditioner + lduVectorPreconditioner; + defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix); + + typedef LduMatrix::smoother + lduVectorSmoother; + defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix); + + typedef LduMatrix::solver + lduVectorSolver; + defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix); + + typedef TPCG TPCGVector; + defineNamedTemplateTypeNameAndDebug(TPCGVector, 0); + + LduMatrix::solver:: + addsymMatrixConstructorToTable + addTPCGSymMatrixConstructorToTable_; + + typedef TPBiCG TPBiCGVector; + defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0); + + LduMatrix::solver:: + addasymMatrixConstructorToTable + addTPBiCGSymMatrixConstructorToTable_; + + typedef NoPreconditioner NoPreconditionerVector; + defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0); + + LduMatrix::preconditioner:: + addsymMatrixConstructorToTable + addNoPreconditionerSymMatrixConstructorToTable_; + + LduMatrix::preconditioner:: + addasymMatrixConstructorToTable + addNoPreconditionerAsymMatrixConstructorToTable_; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Main program: + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + + volVectorField psi + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + lduVectorMatrix testMatrix(mesh); + testMatrix.diag() = 2*pTraits::one; + testMatrix.source() = pTraits::one; + testMatrix.upper() = 0.1; + testMatrix.lower() = -0.1; + + Info<< testMatrix << endl; + + FieldField boundaryCoeffs(0); + FieldField internalCoeffs(0); + + autoPtr testMatrixSolver = + lduVectorMatrix::solver::New + ( + psi.name(), + testMatrix, + boundaryCoeffs, + internalCoeffs, + psi.boundaryField().interfaces(), + IStringStream + ( + "PBiCG" + "{" + " preconditioner none;" + " tolerance (1e-05 1e-05 1e-05);" + " relTol (0 0 0);" + "}" + )() + ); + + lduVectorMatrix::solverPerformance solverPerf = + testMatrixSolver->solve(psi); + + solverPerf.print(); + + Info<< psi << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/LduMatrix/LduMatrixTest2.C b/applications/test/LduMatrix/LduMatrixTest2.C new file mode 100644 index 0000000000..f8694eb481 --- /dev/null +++ b/applications/test/LduMatrix/LduMatrixTest2.C @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 . + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "fvMesh.H" +#include "volFields.H" +#include "LduMatrix.H" +#include "tensorField.H" +#include "TPCG.H" +#include "TPBiCG.H" +#include "NoPreconditioner.H" + +using namespace Foam; + +typedef Foam::LduMatrix + lduVectorMatrix; +defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0); + +typedef Foam::DiagonalSolver + lduVectorDiagonalSolver; +defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0); + +template<> +const vector lduVectorMatrix::great_(1e15, 1e15, 1e15); + +template<> +const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15); + +namespace Foam +{ + typedef LduMatrix::preconditioner + lduVectorPreconditioner; + defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix); + + typedef LduMatrix::smoother + lduVectorSmoother; + defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix); + + typedef LduMatrix::solver + lduVectorSolver; + defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix); + defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix); + + typedef TPCG TPCGVector; + defineNamedTemplateTypeNameAndDebug(TPCGVector, 0); + + LduMatrix::solver:: + addsymMatrixConstructorToTable + addTPCGSymMatrixConstructorToTable_; + + typedef TPBiCG TPBiCGVector; + defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0); + + LduMatrix::solver:: + addasymMatrixConstructorToTable + addTPBiCGSymMatrixConstructorToTable_; + + typedef NoPreconditioner NoPreconditionerVector; + defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0); + + LduMatrix::preconditioner:: + addsymMatrixConstructorToTable + addNoPreconditionerSymMatrixConstructorToTable_; + + LduMatrix::preconditioner:: + addasymMatrixConstructorToTable + addNoPreconditionerAsymMatrixConstructorToTable_; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Main program: + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + + volVectorField psi + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + lduVectorMatrix testMatrix(mesh); + testMatrix.diag() = 2*I; + testMatrix.source() = pTraits::one; + testMatrix.upper() = 0.1; + testMatrix.lower() = -0.1; + + Info<< testMatrix << endl; + + FieldField boundaryCoeffs(0); + FieldField internalCoeffs(0); + + autoPtr testMatrixSolver = + lduVectorMatrix::solver::New + ( + psi.name(), + testMatrix, + //boundaryCoeffs, + //internalCoeffs, + //psi.boundaryField().interfaces(), + IStringStream + ( + "PBiCG" + "{" + " preconditioner none;" + " tolerance (1e-05 1e-05 1e-05);" + " relTol (0 0 0);" + "}" + )() + ); + + lduVectorMatrix::solverPerformance solverPerf = + testMatrixSolver->solve(psi); + + solverPerf.print(); + + Info<< psi << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/LduMatrix/LduMatrixTest3.C b/applications/test/LduMatrix/LduMatrixTest3.C new file mode 100644 index 0000000000..596559e6a7 --- /dev/null +++ b/applications/test/LduMatrix/LduMatrixTest3.C @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 + icoFoam + +Description + Transient solver for incompressible, laminar flow of Newtonian fluids. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "LduMatrix.H" +#include "diagTensorField.H" + +typedef LduMatrix lduVectorMatrix; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +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" + + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + - fvm::laplacian(nu, U) + ); + + fvVectorMatrix UEqnp(UEqn == -fvc::grad(p)); + + lduVectorMatrix U3Eqnp(mesh); + U3Eqnp.diag() = UEqnp.diag(); + U3Eqnp.upper() = UEqnp.upper(); + U3Eqnp.lower() = UEqnp.lower(); + U3Eqnp.source() = UEqnp.source(); + + UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0); + UEqnp.addBoundarySource(U3Eqnp.source(), false); + + autoPtr U3EqnpSolver = + lduVectorMatrix::solver::New + ( + U.name(), + U3Eqnp, + dictionary + ( + IStringStream + ( + "{" + " solver PBiCG;" + " preconditioner DILU;" + " tolerance (1e-5 1e-5 1);" + " relTol (0 0 0);" + "}" + )() + ) + ); + + U3EqnpSolver->solve(U).print(Info); + + // --- PISO loop + + for (int corr=0; corr. + +Application + pisoFoam + +Description + Transient solver for incompressible flow. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulenceModel.H" + +#include "LduMatrix.H" +#include "diagTensorField.H" + +typedef LduMatrix lduVectorMatrix; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +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" + + // Pressure-velocity PISO corrector + { + // Momentum predictor + + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + + turbulence->divDevReff(U) + ); + + //UEqn.relax(); + + fvVectorMatrix UEqnp(UEqn == -fvc::grad(p)); + + lduVectorMatrix U3Eqnp(mesh); + U3Eqnp.diag() = UEqnp.diag(); + U3Eqnp.upper() = UEqnp.upper(); + U3Eqnp.lower() = UEqnp.lower(); + U3Eqnp.source() = UEqnp.source(); + + UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0); + UEqnp.addBoundarySource(U3Eqnp.source(), false); + + U3Eqnp.interfaces() = U.boundaryField().interfaces(); + U3Eqnp.interfacesUpper() = UEqnp.boundaryCoeffs().component(0); + U3Eqnp.interfacesLower() = UEqnp.internalCoeffs().component(0); + + autoPtr U3EqnpSolver = + lduVectorMatrix::solver::New + ( + U.name(), + U3Eqnp, + dictionary + ( + IStringStream + ( + "{" + " /*solver SmoothSolver;*/" + " smoother GaussSeidel;" + " solver PBiCCCG;" + " preconditioner none;" + " tolerance (1e-7 1e-7 1);" + " relTol (0 0 0);" + "}" + )() + ) + ); + + //for (int i=0; i<3; i++) + { + U3EqnpSolver->solve(U).print(Info); + U.correctBoundaryConditions(); + } + //solve(UEqnp); + + // --- PISO loop + + for (int corr=0; corrcorrect(); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/PisoFoam/createFields.H b/applications/test/PisoFoam/createFields.H new file mode 100644 index 0000000000..7cae304f6f --- /dev/null +++ b/applications/test/PisoFoam/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("PISO"), pRefCell, pRefValue); + + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); diff --git a/applications/test/extendedStencil/Test-ExtendedStencil2.C b/applications/test/extendedStencil/Test-ExtendedStencil2.C new file mode 100644 index 0000000000..eda8f98679 --- /dev/null +++ b/applications/test/extendedStencil/Test-ExtendedStencil2.C @@ -0,0 +1,206 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + testExtendedStencil2 + +Description + Test app for determining extended cell-to-cell stencil. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "fvMesh.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "Time.H" +#include "OFstream.H" +#include "meshTools.H" + +#include "CFCCellToCellStencil.H" + + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void writeStencilOBJ +( + const fileName& fName, + const point& fc, + const List& stencilCc +) +{ + OFstream str(fName); + label vertI = 0; + + meshTools::writeOBJ(str, fc); + vertI++; + + forAll(stencilCc, i) + { + meshTools::writeOBJ(str, stencilCc[i]); + vertI++; + str << "l 1 " << vertI << nl; + } +} + + +// Stats +void writeStencilStats(const labelListList& stencil) +{ + label sumSize = 0; + label nSum = 0; + label minSize = labelMax; + label maxSize = labelMin; + + forAll(stencil, i) + { + const labelList& sCells = stencil[i]; + + if (sCells.size() > 0) + { + sumSize += sCells.size(); + nSum++; + minSize = min(minSize, sCells.size()); + maxSize = max(maxSize, sCells.size()); + } + } + reduce(sumSize, sumOp