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/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/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/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/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/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/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/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/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/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/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/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