diff --git a/applications/solvers/combustion/fireFoam/createFieldRefs.H b/applications/solvers/combustion/fireFoam/createFieldRefs.H index 6df6f9df2f..a911c586d9 100644 --- a/applications/solvers/combustion/fireFoam/createFieldRefs.H +++ b/applications/solvers/combustion/fireFoam/createFieldRefs.H @@ -1,4 +1,4 @@ const volScalarField& psi = thermo.psi(); const volScalarField& T = thermo.T(); -filmModelType& surfaceFilm = tsurfaceFilm(); +regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm(); const label inertIndex(composition.species()[inertSpecie]); diff --git a/applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H b/applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H index ffdbcbf6a9..81995c09a5 100644 --- a/applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H +++ b/applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H @@ -1,5 +1,6 @@ Info<< "\nConstructing surface film model" << endl; -typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType; - -autoPtr tsurfaceFilm(filmModelType::New(mesh, g)); +autoPtr tsurfaceFilm +( + regionModels::surfaceFilmModel::New(mesh, g) +); diff --git a/applications/solvers/electromagnetics/mhdFoam/createPhiB.H b/applications/solvers/electromagnetics/mhdFoam/createPhiB.H index 21076a51a5..02b48f8d28 100644 --- a/applications/solvers/electromagnetics/mhdFoam/createPhiB.H +++ b/applications/solvers/electromagnetics/mhdFoam/createPhiB.H @@ -9,7 +9,7 @@ IOobject phiBHeader surfaceScalarField* phiBPtr = nullptr; - if (phiBHeader.typeHeaderOk(true)) +if (phiBHeader.typeHeaderOk(true)) { Info<< "Reading face flux "; diff --git a/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H b/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H index c760349a6e..e7da9c6af6 100644 --- a/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H +++ b/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H @@ -11,7 +11,7 @@ IOobject turbulencePropertiesHeader false ); -if (turbulencePropertiesHeader.typeHeaderOk(false)) +if (turbulencePropertiesHeader.typeHeaderOk(true)) { autoPtr turbulence ( diff --git a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C index aec4c64304..ded2af139b 100644 --- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C +++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C @@ -77,6 +77,9 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; + // Store the particle positions + kinematicCloud.storeGlobalPositions(); + mesh.update(); // Calculate absolute flux from the mapped surface velocity diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C index be0905796b..00c66823ea 100644 --- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C +++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,6 +68,8 @@ int main(int argc, char *argv[]) { Info<< "Time = " << runTime.timeName() << nl << endl; + kinematicCloud.storeGlobalPositions(); + mesh.update(); U.correctBoundaryConditions(); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files deleted file mode 100644 index a76d2acc65..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -reactingParcelFilmFoam.C - -EXE = $(FOAM_APPBIN)/reactingParcelFilmFoam diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H deleted file mode 100644 index 0bb6ff7240..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H +++ /dev/null @@ -1,34 +0,0 @@ - MRF.correctBoundaryVelocity(U); - - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) + fvm::div(phi, U) - + MRF.DDt(rho, U) - + turbulence->divDevRhoReff(U) - == - parcels.SU(U) - + fvOptions(rho, U) - ); - - UEqn.relax(); - - fvOptions.constrain(UEqn); - - if (pimple.momentumPredictor()) - { - solve - ( - UEqn - == - fvc::reconstruct - ( - ( - - ghf*fvc::snGrad(rho) - - fvc::snGrad(p_rgh) - )*mesh.magSf() - ) - ); - - fvOptions.correct(U); - K = 0.5*magSqr(U); - } diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H deleted file mode 100644 index c568be12a1..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H +++ /dev/null @@ -1,9 +0,0 @@ -Info<< "\nConstructing reacting cloud" << endl; -basicReactingCloud parcels -( - "reactingCloud1", - rho, - U, - g, - slgThermo -); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H deleted file mode 100644 index 6ecb8c7e44..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H +++ /dev/null @@ -1,141 +0,0 @@ -Info<< "Creating combustion model\n" << endl; - -autoPtr combustion -( - combustionModels::psiCombustionModel::New(mesh) -); - -psiReactionThermo& thermo = combustion->thermo(); -thermo.validate(args.executable(), "h", "e"); - -SLGThermo slgThermo(mesh, thermo); - -basicSpecieMixture& composition = thermo.composition(); -PtrList& Y = composition.Y(); - -const word inertSpecie(thermo.lookup("inertSpecie")); -if (!composition.species().found(inertSpecie)) -{ - FatalIOErrorIn(args.executable().c_str(), thermo) - << "Inert specie " << inertSpecie << " not found in available species " - << composition.species() - << exit(FatalIOError); -} - -Info<< "Creating field rho\n" << endl; -volScalarField rho -( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - thermo.rho() -); - -volScalarField& p = thermo.p(); - -Info<< "\nReading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "compressibleCreatePhi.H" - -Info<< "Creating turbulence model\n" << endl; -autoPtr turbulence -( - compressible::turbulenceModel::New - ( - rho, - U, - phi, - thermo - ) -); - -// Set the turbulence into the combustion model -combustion->setTurbulence(turbulence()); - -#include "readGravitationalAcceleration.H" -#include "readhRef.H" -#include "gh.H" - - -volScalarField p_rgh -( - IOobject - ( - "p_rgh", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -// Force p_rgh to be consistent with p -p_rgh = p - rho*gh; - -mesh.setFluxRequired(p_rgh.name()); - -multivariateSurfaceInterpolationScheme::fieldTable fields; - -forAll(Y, i) -{ - fields.add(Y[i]); -} -fields.add(thermo.he()); - -IOdictionary additionalControlsDict -( - IOobject - ( - "additionalControls", - runTime.constant(), - mesh, - IOobject::MUST_READ_IF_MODIFIED, - IOobject::NO_WRITE - ) -); - -Switch solvePrimaryRegion -( - additionalControlsDict.lookup("solvePrimaryRegion") -); - -volScalarField Qdot -( - IOobject - ( - "Qdot", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0) -); - -#include "createDpdt.H" - -#include "createK.H" - -#include "createMRF.H" -#include "createClouds.H" -#include "createRadiationModel.H" -#include "createSurfaceFilmModel.H" diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H deleted file mode 100644 index ffdbcbf6a9..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H +++ /dev/null @@ -1,5 +0,0 @@ -Info<< "\nConstructing surface film model" << endl; - -typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType; - -autoPtr tsurfaceFilm(filmModelType::New(mesh, g)); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H deleted file mode 100644 index 7b5249d57e..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H +++ /dev/null @@ -1,59 +0,0 @@ -rho = thermo.rho(); - -volScalarField rAU(1.0/UEqn.A()); -surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); -volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); - -surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - -surfaceScalarField phiHbyA -( - "phiHbyA", - ( - fvc::flux(rho*HbyA) - + rhorAUf*fvc::ddtCorr(rho, U, phi) - ) - + phig -); - -MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - -// Update the pressure BCs to ensure flux consistency -constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); - -while (pimple.correctNonOrthogonal()) -{ - fvScalarMatrix p_rghEqn - ( - fvc::ddt(psi, rho)*gh - + fvc::div(phiHbyA) - + fvm::ddt(psi, p_rgh) - - fvm::laplacian(rhorAUf, p_rgh) - == - parcels.Srho() - + surfaceFilm.Srho() - + fvOptions(psi, p_rgh, rho.name()) - ); - - p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); - - if (pimple.finalNonOrthogonalIter()) - { - phi = phiHbyA + p_rghEqn.flux(); - U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); - U.correctBoundaryConditions(); - fvOptions.correct(U); - } -} - -p = p_rgh + rho*gh; - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -K = 0.5*magSqr(U); - -if (thermo.dpdt()) -{ - dpdt = fvc::ddt(p); -} diff --git a/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H index 42185b8afb..4d112da030 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H @@ -19,6 +19,7 @@ == rho*(U&g) + parcels.Sh(he) + + surfaceFilm.Sh() + radiation->Sh(thermo, he) + Qdot + fvOptions(rho, he) @@ -35,6 +36,6 @@ thermo.correct(); radiation->correct(); - Info<< "T gas min/max " << min(T).value() << ", " + Info<< "T gas min/max = " << min(T).value() << ", " << max(T).value() << endl; } diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options index 8038347006..ca89bd0c01 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options +++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options @@ -1,12 +1,11 @@ EXE_INC = \ -I. \ + -I../reactingParcelFoam \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I${LIB_SRC}/sampling/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ - -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ - -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \ -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/transportModels/compressible/lnInclude \ @@ -17,33 +16,33 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ - -I$(LIB_SRC)/ODE/lnInclude \ -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ -I$(LIB_SRC)/combustionModels/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude \ -I$(FOAM_SOLVERS)/combustion/reactingFoam - EXE_LIBS = \ -lfiniteVolume \ + -lfvOptions \ + -lsampling \ -lmeshTools \ -lturbulenceModels \ -lcompressibleTurbulenceModels \ - -llagrangian \ - -llagrangianIntermediate \ - -llagrangianTurbulence \ -lspecie \ -lcompressibleTransportModels \ -lfluidThermophysicalModels \ - -lthermophysicalProperties \ -lreactionThermophysicalModels \ -lSLGThermo \ -lchemistryModel \ - -lradiationModels \ - -lODE \ -lregionModels \ + -lradiationModels \ -lsurfaceFilmModels \ - -lcombustionModels \ - -lfvOptions \ - -lsampling + -lsurfaceFilmDerivedFvPatchFields \ + -llagrangian \ + -llagrangianIntermediate \ + -llagrangianTurbulence \ + -lODE \ + -lcombustionModels diff --git a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H index bb7ad20efc..0bb6ff7240 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H @@ -6,8 +6,7 @@ + MRF.DDt(rho, U) + turbulence->divDevRhoReff(U) == - rho()*g - + parcels.SU(U) + parcels.SU(U) + fvOptions(rho, U) ); @@ -17,7 +16,18 @@ if (pimple.momentumPredictor()) { - solve(UEqn == -fvc::grad(p)); + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ) + ); fvOptions.correct(U); K = 0.5*magSqr(U); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H index 60a27fec85..eb24899013 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H @@ -9,6 +9,7 @@ tmp> mvConvection ) ); + { combustion->correct(); Qdot = combustion->Qdot(); @@ -24,11 +25,12 @@ tmp> mvConvection ( fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->muEff(), Yi) - == + - fvm::laplacian(turbulence->alphaEff(), Yi) + == parcels.SYi(i, Yi) - + combustion->R(Yi) + fvOptions(rho, Yi) + + combustion->R(Yi) + + surfaceFilm.Srho(i) ); YEqn.relax(); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H b/applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H index e7f76d8456..bbbc50c122 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H @@ -1,3 +1,5 @@ +const label inertIndex(composition.species()[inertSpecie]); + const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); -const label inertIndex(composition.species()[inertSpecie]); +regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm(); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H index ce7e5441e0..b6673dae16 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H @@ -1,7 +1,5 @@ #include "createRDeltaT.H" -#include "readGravitationalAcceleration.H" - Info<< "Creating combustion model\n" << endl; autoPtr combustion @@ -26,8 +24,7 @@ if (!composition.species().found(inertSpecie)) << exit(FatalIOError); } -volScalarField& p = thermo.p(); - +Info<< "Creating field rho\n" << endl; volScalarField rho ( IOobject @@ -41,6 +38,8 @@ volScalarField rho thermo.rho() ); +volScalarField& p = thermo.p(); + Info<< "\nReading field U\n" << endl; volVectorField U ( @@ -57,30 +56,6 @@ volVectorField U #include "compressibleCreatePhi.H" -mesh.setFluxRequired(p.name()); - -dimensionedScalar rhoMax -( - dimensionedScalar::lookupOrDefault - ( - "rhoMax", - pimple.dict(), - dimDensity, - GREAT - ) -); - -dimensionedScalar rhoMin -( - dimensionedScalar::lookupOrDefault - ( - "rhoMin", - pimple.dict(), - dimDensity, - 0 - ) -); - Info<< "Creating turbulence model\n" << endl; autoPtr turbulence ( @@ -96,6 +71,31 @@ autoPtr turbulence // Set the turbulence into the combustion model combustion->setTurbulence(turbulence()); +#include "readGravitationalAcceleration.H" +#include "readhRef.H" +#include "gh.H" + + +volScalarField p_rgh +( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +// Force p_rgh to be consistent with p +p_rgh = p - rho*gh; + +pressureControl pressureControl(p, rho, pimple.dict(), false); + +mesh.setFluxRequired(p_rgh.name()); + Info<< "Creating multi-variate interpolation scheme\n" << endl; multivariateSurfaceInterpolationScheme::fieldTable fields; @@ -105,6 +105,11 @@ forAll(Y, i) } fields.add(thermo.he()); +Switch solvePrimaryRegion +( + pimple.dict().lookupOrDefault("solvePrimaryRegion", true) +); + volScalarField Qdot ( IOobject @@ -126,3 +131,4 @@ volScalarField Qdot #include "createMRF.H" #include "createRadiationModel.H" #include "createClouds.H" +#include "createSurfaceFilmModel.H" diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createSurfaceFilmModel.H b/applications/solvers/lagrangian/reactingParcelFoam/createSurfaceFilmModel.H new file mode 100644 index 0000000000..81995c09a5 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/createSurfaceFilmModel.H @@ -0,0 +1,6 @@ +Info<< "\nConstructing surface film model" << endl; + +autoPtr tsurfaceFilm +( + regionModels::surfaceFilmModel::New(mesh, g) +); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H index f2cb294718..588273cf54 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H @@ -1,4 +1,7 @@ -rho = thermo.rho(); +if (!pimple.SIMPLErho()) +{ + rho = thermo.rho(); +} // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution @@ -7,6 +10,9 @@ const volScalarField psip0(psi*p); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + surfaceScalarField phiHbyA ( "phiHbyA", @@ -14,58 +20,68 @@ surfaceScalarField phiHbyA fvc::flux(rho*HbyA) + rhorAUf*fvc::ddtCorr(rho, U, phi) ) + + phig ); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); // Update the pressure BCs to ensure flux consistency -constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); +constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); -fvScalarMatrix pDDtEqn +fvScalarMatrix p_rghDDtEqn ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phiHbyA) == parcels.Srho() - + fvOptions(psi, p, rho.name()) + + surfaceFilm.Srho() + + fvOptions(psi, p_rgh, rho.name()) ); while (pimple.correctNonOrthogonal()) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - pDDtEqn - - fvm::laplacian(rhorAUf, p) + p_rghDDtEqn + - fvm::laplacian(rhorAUf, p_rgh) ); - pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { - phi = phiHbyA + pEqn.flux(); + phi = phiHbyA + p_rghEqn.flux(); + + // Explicitly relax pressure for momentum corrector + p_rgh.relax(); + + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); + U.correctBoundaryConditions(); + fvOptions.correct(U); + K = 0.5*magSqr(U); } } -p.relax(); +p = p_rgh + rho*gh; // Thermodynamic density update thermo.correctRho(psi*p - psip0); -#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent +#include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U = HbyA - rAU*fvc::grad(p); -U.correctBoundaryConditions(); -fvOptions.correct(U); -K = 0.5*magSqr(U); +if (pressureControl.limit(p)) +{ + p.correctBoundaryConditions(); + rho = thermo.rho(); + p_rgh = p - rho*gh; +} +else if (pimple.SIMPLErho()) +{ + rho = thermo.rho(); +} if (thermo.dpdt()) { dpdt = fvc::ddt(p); } - -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); - -Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl; diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C index 93768896b3..497aa835a5 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,18 +29,20 @@ Group Description Transient solver for compressible, turbulent flow with a reacting, - multiphase particle cloud, and optional sources/constraints. + multiphase particle cloud, and surface film modelling. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "turbulentFluidThermoModel.H" #include "basicReactingMultiphaseCloud.H" +#include "surfaceFilmModel.H" #include "rhoCombustionModel.H" #include "radiationModel.H" -#include "fvOptions.H" #include "SLGThermo.H" +#include "fvOptions.H" #include "pimpleControl.H" +#include "pressureControl.H" #include "localEulerDdtScheme.H" #include "fvcSmooth.H" @@ -76,10 +78,14 @@ int main(int argc, char *argv[]) { #include "readTimeControls.H" - if (!LTS) + if (LTS) + { + #include "setRDeltaT.H" + } + else { #include "compressibleCourantNo.H" - #include "setDeltaT.H" + #include "setMultiRegionDeltaT.H" } runTime++; @@ -87,34 +93,36 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; parcels.evolve(); + surfaceFilm.evolve(); - if (LTS) + if (solvePrimaryRegion) { - #include "setRDeltaT.H" - } - - #include "rhoEqn.H" - - // --- Pressure-velocity PIMPLE corrector loop - while (pimple.loop()) - { - #include "UEqn.H" - #include "YEqn.H" - #include "EEqn.H" - - // --- Pressure corrector loop - while (pimple.correct()) + if (pimple.nCorrPIMPLE() <= 1) { - #include "pEqn.H" + #include "rhoEqn.H" } - if (pimple.turbCorr()) + // --- PIMPLE loop + while (pimple.loop()) { - turbulence->correct(); - } - } + #include "UEqn.H" + #include "YEqn.H" + #include "EEqn.H" - rho = thermo.rho(); + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + rho = thermo.rho(); + } runTime.write(); @@ -123,7 +131,7 @@ int main(int argc, char *argv[]) << nl << endl; } - Info<< "End\n" << endl; + Info<< "End" << endl; return 0; } diff --git a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H index b69701e743..baf1cb4ec4 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,15 +36,13 @@ Description + fvc::div(phi) == parcels.Srho(rho) + + surfaceFilm.Srho() + fvOptions(rho) ); rhoEqn.solve(); fvOptions.correct(rho); - - Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value() - << endl; } // ************************************************************************* // diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H similarity index 79% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H rename to applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H index 76d08256eb..a61829a6be 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H @@ -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-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -33,25 +33,21 @@ Description if (adjustTimeStep) { - if (CoNum == -GREAT) - { - CoNum = SMALL; - } - - const scalar TFactorFluid = maxCo/(CoNum + SMALL); - const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL); - - const scalar dt0 = runTime.deltaTValue(); + const scalar maxDeltaTFact = + min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL)); + const scalar deltaTFact = + min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); runTime.setDeltaT ( min ( - dt0*min(min(TFactorFluid, TFactorFilm), 1.2), + deltaTFact*runTime.deltaTValue(), maxDeltaT ) ); + + Info<< "deltaT = " << runTime.deltaTValue() << endl; } - // ************************************************************************* // diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C index bbec4a41ab..a5092143a9 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,9 +24,6 @@ License Application simpleReactingParcelFoam -Group - grpLagrangianSolvers - Description Steady state solver for compressible, turbulent flow with reacting, multiphase particle clouds and optional sources/constraints. @@ -38,6 +35,7 @@ Description #include "basicReactingMultiphaseCloud.H" #include "rhoCombustionModel.H" #include "radiationModel.H" +#include "IOporosityModelList.H" #include "fvOptions.H" #include "SLGThermo.H" #include "simpleControl.H" diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H b/applications/solvers/lagrangian/sprayFoam/EEqn.H similarity index 90% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H rename to applications/solvers/lagrangian/sprayFoam/EEqn.H index 4d112da030..42185b8afb 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/EEqn.H @@ -19,7 +19,6 @@ == rho*(U&g) + parcels.Sh(he) - + surfaceFilm.Sh() + radiation->Sh(thermo, he) + Qdot + fvOptions(rho, he) @@ -36,6 +35,6 @@ thermo.correct(); radiation->correct(); - Info<< "T gas min/max = " << min(T).value() << ", " + Info<< "T gas min/max " << min(T).value() << ", " << max(T).value() << endl; } diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/sprayFoam/YEqn.H similarity index 89% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H rename to applications/solvers/lagrangian/sprayFoam/YEqn.H index eb24899013..60a27fec85 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/YEqn.H @@ -9,7 +9,6 @@ tmp> mvConvection ) ); - { combustion->correct(); Qdot = combustion->Qdot(); @@ -25,12 +24,11 @@ tmp> mvConvection ( fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->alphaEff(), Yi) - == + - fvm::laplacian(turbulence->muEff(), Yi) + == parcels.SYi(i, Yi) - + fvOptions(rho, Yi) + combustion->R(Yi) - + surfaceFilm.Srho(i) + + fvOptions(rho, Yi) ); YEqn.relax(); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFieldRefs.H b/applications/solvers/lagrangian/sprayFoam/createFieldRefs.H similarity index 75% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/createFieldRefs.H rename to applications/solvers/lagrangian/sprayFoam/createFieldRefs.H index 14914688ca..e7f76d8456 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFieldRefs.H +++ b/applications/solvers/lagrangian/sprayFoam/createFieldRefs.H @@ -1,5 +1,3 @@ -const label inertIndex(composition.species()[inertSpecie]); - const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); -filmModelType& surfaceFilm = tsurfaceFilm(); +const label inertIndex(composition.species()[inertSpecie]); diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C index 16d2d524c5..40acccae9e 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C +++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,6 +91,9 @@ int main(int argc, char *argv[]) // Store momentum to set rhoUf for introduced faces. volVectorField rhoU("rhoU", rho*U); + // Store the particle positions + parcels.storeGlobalPositions(); + // Do any mesh changes mesh.update(); diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files new file mode 100644 index 0000000000..430b1df24f --- /dev/null +++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files @@ -0,0 +1,3 @@ +uncoupledKinematicParcelDyMFoam.C + +EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelDyMFoam diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options similarity index 59% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options rename to applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options index eab72a2389..2d385f9853 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options +++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options @@ -1,46 +1,36 @@ EXE_INC = \ - -I. \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I${LIB_SRC}/sampling/lnInclude \ - -I${LIB_SRC}/meshTools/lnInclude \ - -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ - -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ - -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ + -I.. \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/transportModels/compressible/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ - -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ - -I$(LIB_SRC)/ODE/lnInclude \ - -I$(LIB_SRC)/combustionModels/lnInclude \ - -I$(FOAM_SOLVERS)/combustion/reactingFoam + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude EXE_LIBS = \ - -lfiniteVolume \ - -lfvOptions \ - -lsampling \ - -lmeshTools \ - -lturbulenceModels \ - -lcompressibleTurbulenceModels \ - -lspecie \ - -lcompressibleTransportModels \ - -lfluidThermophysicalModels \ - -lreactionThermophysicalModels \ - -lSLGThermo \ - -lchemistryModel \ - -lregionModels \ - -lradiationModels \ - -lsurfaceFilmModels \ - -lsurfaceFilmDerivedFvPatchFields \ + -llagrangian \ -llagrangianIntermediate \ -llagrangianTurbulence \ - -lODE \ - -lcombustionModels + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lradiationModels \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lregionModels \ + -lsurfaceFilmModels \ + -ldynamicMesh \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C similarity index 54% rename from applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C rename to applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C index 35348ac933..18f8e63869 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C +++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,91 +22,57 @@ License along with OpenFOAM. If not, see . Application - reactingParcelFilmFoam - -Group - grpLagrangianSolvers + uncoupledKinematicParcelDyMFoam Description - Transient solver for compressible, turbulent flow with a reacting, - multiphase particle cloud, and surface film modelling. + Transient solver for the passive transport of a particle cloud. + + Uses a pre- calculated velocity field to evolve the cloud. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "psiThermo.H" #include "turbulentFluidThermoModel.H" -#include "basicReactingCloud.H" -#include "surfaceFilmModel.H" -#include "psiCombustionModel.H" -#include "radiationModel.H" -#include "SLGThermo.H" -#include "fvOptions.H" -#include "pimpleControl.H" +#include "basicKinematicCloud.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { + argList::addOption + ( + "cloudName", + "name", + "specify alternative cloud name. default is 'kinematicCloud'" + ); + + #define NO_CONTROL #include "postProcess.H" #include "setRootCase.H" #include "createTime.H" - #include "createMesh.H" - #include "createControl.H" - #include "createTimeControls.H" + #include "createDynamicFvMesh.H" #include "createFields.H" - #include "createFieldRefs.H" - #include "createFvOptions.H" #include "compressibleCourantNo.H" - #include "setInitialDeltaT.H" - #include "initContinuityErrs.H" - - turbulence->validate(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; - while (runTime.run()) + while (runTime.loop()) { - #include "readTimeControls.H" - #include "compressibleCourantNo.H" - #include "setMultiRegionDeltaT.H" - #include "setDeltaT.H" - - runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; - parcels.evolve(); + kinematicCloud.storeGlobalPositions(); - surfaceFilm.evolve(); + mesh.update(); - if (solvePrimaryRegion) - { - #include "rhoEqn.H" + U.correctBoundaryConditions(); - // --- PIMPLE loop - while (pimple.loop()) - { - #include "UEqn.H" - #include "YEqn.H" - #include "EEqn.H" - - // --- Pressure corrector loop - while (pimple.correct()) - { - #include "pEqn.H" - } - - if (pimple.turbCorr()) - { - turbulence->correct(); - } - } - - rho = thermo.rho(); - } + Info<< "Evolving " << kinematicCloud.name() << endl; + kinematicCloud.evolve(); runTime.write(); @@ -115,7 +81,7 @@ int main(int argc, char *argv[]) << nl << endl; } - Info<< "End" << endl; + Info<< "End\n" << endl; return 0; } diff --git a/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H index 3ce881bd06..40ca20533c 100644 --- a/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H +++ b/applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H @@ -118,27 +118,27 @@ << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; - tmp talphaPhiUD(alpha1Eqn.flux()); - alphaPhi = talphaPhiUD(); + tmp talphaPhi1UD(alpha1Eqn.flux()); + alphaPhi10 = talphaPhi1UD(); - if (alphaApplyPrevCorr && talphaPhiCorr0.valid()) + if (alphaApplyPrevCorr && talphaPhi1Corr0.valid()) { Info<< "Applying the previous iteration compression flux" << endl; MULES::correct ( alphac, alpha1, - alphaPhi, - talphaPhiCorr0.ref(), + alphaPhi10, + talphaPhi1Corr0.ref(), zeroField(), zeroField(), 1, 0 ); - alphaPhi += talphaPhiCorr0(); + alphaPhi10 += talphaPhi1Corr0(); } // Cache the upwind-flux - talphaPhiCorr0 = talphaPhiUD; + talphaPhi1Corr0 = talphaPhi1UD; alpha2 = 1.0 - alpha1; @@ -152,7 +152,7 @@ surfaceScalarField phir(phic*mixture.nHatf()); - tmp talphaPhiUn + tmp talphaPhi1Un ( fvc::flux ( @@ -170,15 +170,15 @@ if (MULESCorr) { - tmp talphaPhiCorr(talphaPhiUn() - alphaPhi); + tmp talphaPhi1Corr(talphaPhi1Un() - alphaPhi10); volScalarField alpha10("alpha10", alpha1); MULES::correct ( alphac, alpha1, - talphaPhiUn(), - talphaPhiCorr.ref(), + talphaPhi1Un(), + talphaPhi1Corr.ref(), Sp, (-Sp*alpha1)(), 1, @@ -188,24 +188,24 @@ // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { - alphaPhi += talphaPhiCorr(); + alphaPhi10 += talphaPhi1Corr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; - alphaPhi += 0.5*talphaPhiCorr(); + alphaPhi10 += 0.5*talphaPhi1Corr(); } } else { - alphaPhi = talphaPhiUn; + alphaPhi10 = talphaPhi1Un; MULES::explicitSolve ( alphac, alpha1, phiCN, - alphaPhi, + alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), 1, @@ -220,12 +220,12 @@ if (alphaApplyPrevCorr && MULESCorr) { - talphaPhiCorr0 = alphaPhi - talphaPhiCorr0; - talphaPhiCorr0.ref().rename("alphaPhiCorr0"); + talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0; + talphaPhi1Corr0.ref().rename("alphaPhi1Corr0"); } else { - talphaPhiCorr0.clear(); + talphaPhi1Corr0.clear(); } if @@ -235,19 +235,20 @@ ) { #include "rhofs.H" - rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f; + rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f; } else { if (ocCoeff > 0) { // Calculate the end-of-time-step alpha flux - alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff; + alphaPhi10 = + (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff; } // Calculate the end-of-time-step mass flux #include "rhofs.H" - rhoPhi = alphaPhi*(rho1f - rho2f) + alphaPhic*rho2f; + rhoPhi = alphaPhi10*(rho1f - rho2f) + alphaPhic*rho2f; } Info<< "Phase-1 volume fraction = " diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H index abae00fb38..ad9ddcb5e5 100644 --- a/applications/solvers/multiphase/VoF/alphaEqn.H +++ b/applications/solvers/multiphase/VoF/alphaEqn.H @@ -65,6 +65,14 @@ phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } + // Add the optional shear compression contribution + if (scAlpha > 0) + { + phic += + scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(U)))); + } + + surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef(); @@ -105,6 +113,8 @@ phiCN, upwind(mesh, phiCN) ).fvmDiv(phiCN, alpha1) + // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh) + // + fvc::div(phiCN), alpha1) == Su + fvm::Sp(Sp + divU, alpha1) ); @@ -117,19 +127,19 @@ << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; - tmp talphaPhiUD(alpha1Eqn.flux()); - alphaPhi = talphaPhiUD(); + tmp talphaPhi1UD(alpha1Eqn.flux()); + alphaPhi10 = talphaPhi1UD(); - if (alphaApplyPrevCorr && talphaPhiCorr0.valid()) + if (alphaApplyPrevCorr && talphaPhi1Corr0.valid()) { Info<< "Applying the previous iteration compression flux" << endl; - MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0); + MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0); - alphaPhi += talphaPhiCorr0(); + alphaPhi10 += talphaPhi1Corr0(); } // Cache the upwind-flux - talphaPhiCorr0 = talphaPhiUD; + talphaPhi1Corr0 = talphaPhi1UD; alpha2 = 1.0 - alpha1; @@ -143,7 +153,7 @@ surfaceScalarField phir(phic*mixture.nHatf()); - tmp talphaPhiUn + tmp talphaPhi1Un ( fvc::flux ( @@ -161,15 +171,15 @@ if (MULESCorr) { - tmp talphaPhiCorr(talphaPhiUn() - alphaPhi); + tmp talphaPhi1Corr(talphaPhi1Un() - alphaPhi10); volScalarField alpha10("alpha10", alpha1); MULES::correct ( geometricOneField(), alpha1, - talphaPhiUn(), - talphaPhiCorr.ref(), + talphaPhi1Un(), + talphaPhi1Corr.ref(), Sp, (-Sp*alpha1)(), 1, @@ -179,24 +189,24 @@ // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { - alphaPhi += talphaPhiCorr(); + alphaPhi10 += talphaPhi1Corr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; - alphaPhi += 0.5*talphaPhiCorr(); + alphaPhi10 += 0.5*talphaPhi1Corr(); } } else { - alphaPhi = talphaPhiUn; + alphaPhi10 = talphaPhi1Un; MULES::explicitSolve ( geometricOneField(), alpha1, phiCN, - alphaPhi, + alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), 1, @@ -211,34 +221,37 @@ if (alphaApplyPrevCorr && MULESCorr) { - talphaPhiCorr0 = alphaPhi - talphaPhiCorr0; - talphaPhiCorr0.ref().rename("alphaPhiCorr0"); + talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0; + talphaPhi1Corr0.ref().rename("alphaPhi1Corr0"); } else { - talphaPhiCorr0.clear(); + talphaPhi1Corr0.clear(); } + #include "rhofs.H" + if ( word(mesh.ddtScheme("ddt(rho,U)")) == fv::EulerDdtScheme::typeName + || word(mesh.ddtScheme("ddt(rho,U)")) + == fv::localEulerDdtScheme::typeName ) { - #include "rhofs.H" - rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f; + rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f; } else { if (ocCoeff > 0) { // Calculate the end-of-time-step alpha flux - alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff; + alphaPhi10 = + (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff; } // Calculate the end-of-time-step mass flux - #include "rhofs.H" - rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f; + rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f; } Info<< "Phase-1 volume fraction = " diff --git a/applications/solvers/multiphase/VoF/createAlphaFluxes.H b/applications/solvers/multiphase/VoF/createAlphaFluxes.H index 3a697e38d4..a9004f2f8e 100644 --- a/applications/solvers/multiphase/VoF/createAlphaFluxes.H +++ b/applications/solvers/multiphase/VoF/createAlphaFluxes.H @@ -1,20 +1,21 @@ -IOobject alphaPhiHeader +IOobject alphaPhi10Header ( - "alphaPhi", + "alphaPhi10", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ); -const bool alphaRestart = alphaPhiHeader.typeHeaderOk(true); +const bool alphaRestart = + alphaPhi10Header.typeHeaderOk(true); // MULES flux from previous time-step -surfaceScalarField alphaPhi +surfaceScalarField alphaPhi10 ( - alphaPhiHeader, + alphaPhi10Header, phi*fvc::interpolate(alpha1) ); // MULES Correction -tmp talphaPhiCorr0; +tmp talphaPhi1Corr0; diff --git a/applications/solvers/multiphase/VoF/setRDeltaT.H b/applications/solvers/multiphase/VoF/setRDeltaT.H index 697e6588e2..26b237e795 100644 --- a/applications/solvers/multiphase/VoF/setRDeltaT.H +++ b/applications/solvers/multiphase/VoF/setRDeltaT.H @@ -73,8 +73,8 @@ rDeltaT.ref() = max ( rDeltaT(), - pos(alpha1Bar() - alphaSpreadMin) - *pos(alphaSpreadMax - alpha1Bar()) + pos0(alpha1Bar() - alphaSpreadMin) + *pos0(alphaSpreadMax - alpha1Bar()) *fvc::surfaceSum(mag(phi))()() /((2*maxAlphaCo)*mesh.V()) ); diff --git a/applications/solvers/multiphase/compressibleInterFoam/Allwclean b/applications/solvers/multiphase/compressibleInterFoam/Allwclean index 71bff64a72..f06d2f228d 100755 --- a/applications/solvers/multiphase/compressibleInterFoam/Allwclean +++ b/applications/solvers/multiphase/compressibleInterFoam/Allwclean @@ -5,5 +5,6 @@ wclean libso twoPhaseMixtureThermo wclean libso surfaceTensionModels wclean wclean compressibleInterDyMFoam +wclean compressibleInterFilmFoam #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/compressibleInterFoam/Allwmake b/applications/solvers/multiphase/compressibleInterFoam/Allwmake index 29bafc3a61..7ecc9f4238 100755 --- a/applications/solvers/multiphase/compressibleInterFoam/Allwmake +++ b/applications/solvers/multiphase/compressibleInterFoam/Allwmake @@ -9,5 +9,6 @@ wmake $targetType surfaceTensionModels wmake $targetType wmake $targetType compressibleInterDyMFoam +wmake $targetType compressibleInterFilmFoam #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H index af74f0362a..9f9e64d199 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H @@ -1,8 +1,8 @@ { fvScalarMatrix TEqn ( - fvm::ddt(rho, T) - + fvm::div(rhoPhi, T) + fvm::ddt(rho, T) + fvm::div(rhoPhi, T) + - fvm::Sp(contErr, T) - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T) + ( divU*p diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H index b5e66ec33f..d5424cc72f 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H @@ -1,6 +1,7 @@ fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + - fvm::Sp(contErr, U) + MRF.DDt(rho, U) + turbulence->divDevRhoReff(U) == diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H index 81309cd091..12d2bcfbd7 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H +++ b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H @@ -24,14 +24,14 @@ volScalarField::Internal Su forAll(dgdt, celli) { - if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0) + if (dgdt[celli] > 0.0) { - Sp[celli] -= dgdt[celli]*alpha1[celli]; - Su[celli] += dgdt[celli]*alpha1[celli]; + Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); + Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); } - else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0) + else if (dgdt[celli] < 0.0) { - Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]); + Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4); } } diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H index f8d0a15836..006b26317e 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H @@ -1,3 +1,69 @@ +tmp talphaPhi1(alphaPhi10); + +if (nAlphaSubCycles > 1) { - #include "alphaEqnSubCycle.H" + dimensionedScalar totalDeltaT = runTime.deltaT(); + + talphaPhi1 = new surfaceScalarField + ( + IOobject + ( + "alphaPhi1", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("0", alphaPhi10.dimensions(), 0) + ); + + surfaceScalarField rhoPhiSum + ( + IOobject + ( + "rhoPhiSum", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("0", rhoPhi.dimensions(), 0) + ); + + tmp trSubDeltaT; + + if (LTS) + { + trSubDeltaT = + fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles); + } + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqn.H" + talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi10; + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; } +else +{ + #include "alphaEqn.H" +} + +rho == alpha1*rho1 + alpha2*rho2; + +const surfaceScalarField& alphaPhi1 = talphaPhi1(); +surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1); + +volScalarField::Internal contErr +( + ( + fvc::ddt(rho) + fvc::div(rhoPhi) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )() +); diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C index e3d43a8944..0366272631 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C @@ -150,8 +150,6 @@ int main(int argc, char *argv[]) #include "alphaControls.H" #include "compressibleAlphaEqnSubCycle.H" - solve(fvm::ddt(rho) + fvc::div(rhoPhi)); - #include "UEqn.H" #include "TEqn.H" diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H index e84c19453f..b1eace14f2 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H @@ -56,13 +56,29 @@ } else { + #include "rhofs.H" + p_rghEqnComp1 = - fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh)) - + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1); + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)) + ); p_rghEqnComp2 = - fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh)) - + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2); + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)) + ); } // Cache p_rgh prior to solve for density update @@ -78,16 +94,21 @@ solve ( - ( - (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1() - + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2() - ) - + p_rghEqnIncomp, + p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp, mesh.solver(p_rgh.select(pimple.finalInnerIter())) ); if (pimple.finalNonOrthogonalIter()) { + p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); + p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; + + dgdt = + ( + alpha1*(p_rghEqnComp2 & p_rgh) + - alpha2*(p_rghEqnComp1 & p_rgh) + ); + phi = phiHbyA + p_rghEqnIncomp.flux(); U = HbyA @@ -109,15 +130,9 @@ rho = alpha1*rho1 + alpha2*rho2; - p = max(p_rgh + rho*gh, pMin); + // Correct p_rgh for consistency with p and the updated densities p_rgh = p - rho*gh; p_rgh.correctBoundaryConditions(); - dgdt = - ( - pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 - - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 - ); - K = 0.5*magSqr(U); } diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/files b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/files new file mode 100644 index 0000000000..d24eb3216c --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/files @@ -0,0 +1,6 @@ +VoFPatchTransfer/VoFPatchTransfer.C +VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C +VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C +compressibleInterFilmFoam.C + +EXE = $(FOAM_APPBIN)/compressibleInterFilmFoam diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options new file mode 100644 index 0000000000..5827d0b9fd --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options @@ -0,0 +1,44 @@ +EXE_INC = \ + -I. \ + -I.. \ + -I../../VoF \ + -I../twoPhaseMixtureThermo \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude + +EXE_LIBS = \ + -ltwoPhaseMixtureThermo \ + -ltwoPhaseSurfaceTension \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -ltwoPhaseMixture \ + -ltwoPhaseProperties \ + -linterfaceProperties \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lSLGThermo \ + -lsurfaceFilmModels \ + -lsurfaceFilmDerivedFvPatchFields \ + -ldynamicMesh \ + -lmeshTools \ + -ldynamicFvMesh \ + -lfiniteVolume \ + -lfvOptions diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/TEqn.H new file mode 100644 index 0000000000..11475fa4bb --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/TEqn.H @@ -0,0 +1,36 @@ +{ + fvScalarMatrix TEqn + ( + fvm::ddt(rho, T) + fvm::div(rhoPhi, T) + - fvm::Sp(contErr, T) + - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T) + + ( + ( + fvc::div(fvc::absolute(phi, U), p) + + fvc::ddt(rho, K) + fvc::div(rhoPhi, K) + )()() - contErr*K + ) + *( + alpha1/mixture.thermo1().Cv() + + alpha2/mixture.thermo2().Cv() + )()() + == + fvOptions(rho, T) + + pos(Srho) + *( + surfaceFilm.Sh()()/mixture.thermo1().Cp()() + + surfaceFilm.Tref*Srho + ) + ); + + TEqn.relax(); + + fvOptions.constrain(TEqn); + + TEqn.solve(); + + fvOptions.correct(T); + + mixture.correctThermo(); + mixture.correct(); +} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C new file mode 100644 index 0000000000..b3a646e93b --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C @@ -0,0 +1,328 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "VoFPatchTransfer.H" +#include "twoPhaseMixtureThermo.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace regionModels +{ +namespace surfaceFilmModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(VoFPatchTransfer, 0); +addToRunTimeSelectionTable(transferModel, VoFPatchTransfer, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +VoFPatchTransfer::VoFPatchTransfer +( + surfaceFilmRegionModel& film, + const dictionary& dict +) +: + transferModel(type(), film, dict), + deltaFactorToVoF_ + ( + coeffDict_.lookupOrDefault("deltaFactorToVoF", 1.0) + ), + deltaFactorToFilm_ + ( + coeffDict_.lookupOrDefault("deltaFactorToFilm", 0.5) + ), + alphaToVoF_ + ( + coeffDict_.lookupOrDefault("alphaToVoF", 0.5) + ), + alphaToFilm_ + ( + coeffDict_.lookupOrDefault("alphaToFilm", 0.1) + ), + transferRateCoeff_ + ( + coeffDict_.lookupOrDefault("transferRateCoeff", 0.1) + ) +{ + const polyBoundaryMesh& pbm = film.regionMesh().boundaryMesh(); + patchIDs_.setSize + ( + pbm.size() - film.regionMesh().globalData().processorPatches().size() + ); + + if (coeffDict_.found("patches")) + { + const wordReList patchNames(coeffDict_.lookup("patches")); + const labelHashSet patchSet = pbm.patchSet(patchNames); + + Info<< " applying to patches:" << nl; + + label pidi = 0; + forAllConstIter(labelHashSet, patchSet, iter) + { + const label patchi = iter.key(); + patchIDs_[pidi++] = patchi; + Info<< " " << pbm[patchi].name() << endl; + } + patchIDs_.setSize(pidi); + patchTransferredMasses_.setSize(pidi, 0); + } + else + { + Info<< " applying to all patches" << endl; + + forAll(patchIDs_, patchi) + { + patchIDs_[patchi] = patchi; + } + + patchTransferredMasses_.setSize(patchIDs_.size(), 0); + } + + if (!patchIDs_.size()) + { + FatalErrorInFunction + << "No patches selected" + << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +VoFPatchTransfer::~VoFPatchTransfer() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void VoFPatchTransfer::correct +( + scalarField& availableMass, + scalarField& massToTransfer +) +{ + NotImplemented; +} + + +void VoFPatchTransfer::correct +( + scalarField& availableMass, + scalarField& massToTransfer, + scalarField& energyToTransfer +) +{ + // Do not correct if no patches selected + if (!patchIDs_.size()) return; + + const scalarField& delta = film().delta(); + const scalarField& rho = film().rho(); + const scalarField& magSf = film().magSf(); + + const polyBoundaryMesh& pbm = film().regionMesh().boundaryMesh(); + + + const twoPhaseMixtureThermo& thermo + ( + film().primaryMesh().lookupObject + ( + twoPhaseMixtureThermo::dictName + ) + ); + + const volScalarField& heVoF = thermo.thermo1().he(); + const volScalarField& TVoF = thermo.thermo1().T(); + const volScalarField CpVoF(thermo.thermo1().Cp()); + const volScalarField& rhoVoF = thermo.thermo1().rho()(); + const volScalarField& alphaVoF = thermo.alpha1(); + + forAll(patchIDs_, pidi) + { + const label patchi = patchIDs_[pidi]; + label primaryPatchi = -1; + + forAll(film().intCoupledPatchIDs(), i) + { + const label filmPatchi = film().intCoupledPatchIDs()[i]; + + if (filmPatchi == patchi) + { + primaryPatchi = film().primaryPatchIDs()[i]; + } + } + + if (primaryPatchi != -1) + { + scalarField deltaCoeffs + ( + film().primaryMesh().boundary()[primaryPatchi].deltaCoeffs() + ); + film().toRegion(patchi, deltaCoeffs); + + scalarField hp(heVoF.boundaryField()[primaryPatchi]); + film().toRegion(patchi, hp); + + scalarField Tp(TVoF.boundaryField()[primaryPatchi]); + film().toRegion(patchi, Tp); + + scalarField Cpp(CpVoF.boundaryField()[primaryPatchi]); + film().toRegion(patchi, Cpp); + + scalarField rhop(rhoVoF.boundaryField()[primaryPatchi]); + film().toRegion(patchi, rhop); + + scalarField alphap(alphaVoF.boundaryField()[primaryPatchi]); + film().toRegion(patchi, alphap); + + scalarField Vp + ( + film().primaryMesh().boundary()[primaryPatchi] + .patchInternalField(film().primaryMesh().V()) + ); + film().toRegion(patchi, Vp); + + const polyPatch& pp = pbm[patchi]; + const labelList& faceCells = pp.faceCells(); + + // Accumulate the total mass removed from patch + scalar dMassPatch = 0; + + forAll(faceCells, facei) + { + const label celli = faceCells[facei]; + + scalar dMass = 0; + + if + ( + delta[celli] > 2*deltaFactorToVoF_/deltaCoeffs[facei] + || alphap[facei] > alphaToVoF_ + ) + { + dMass = + transferRateCoeff_*delta[celli]*rho[celli]*magSf[celli]; + + massToTransfer[celli] += dMass; + energyToTransfer[celli] += dMass*film().hs()[celli]; + } + + if + ( + alphap[facei] > 0 + && delta[celli] > 0 + && delta[celli] < 2*deltaFactorToFilm_/deltaCoeffs[facei] + && alphap[facei] < alphaToFilm_ + ) + { + dMass = + -transferRateCoeff_*alphap[facei]*rhop[facei]*Vp[facei]; + + massToTransfer[celli] += dMass; + energyToTransfer[celli] += dMass*hp[facei]; + } + + availableMass[celli] -= dMass; + dMassPatch += dMass; + } + + patchTransferredMasses_[pidi] += dMassPatch; + addToTransferredMass(dMassPatch); + } + } + + transferModel::correct(); + + if (writeTime()) + { + scalarField patchTransferredMasses0 + ( + getModelProperty + ( + "patchTransferredMasses", + scalarField(patchTransferredMasses_.size(), 0) + ) + ); + + scalarField patchTransferredMassTotals(patchTransferredMasses_); + Pstream::listCombineGather + ( + patchTransferredMassTotals, + plusEqOp() + ); + patchTransferredMasses0 += patchTransferredMassTotals; + + setModelProperty + ( + "patchTransferredMasses", + patchTransferredMasses0 + ); + + patchTransferredMasses_ = 0; + } +} + + +void VoFPatchTransfer::patchTransferredMassTotals +( + scalarField& patchMasses +) const +{ + // Do not correct if no patches selected + if (!patchIDs_.size()) return; + + scalarField patchTransferredMasses + ( + getModelProperty + ( + "patchTransferredMasses", + scalarField(patchTransferredMasses_.size(), 0) + ) + ); + + scalarField patchTransferredMassTotals(patchTransferredMasses_); + Pstream::listCombineGather(patchTransferredMassTotals, plusEqOp()); + + forAll(patchIDs_, pidi) + { + const label patchi = patchIDs_[pidi]; + patchMasses[patchi] += + patchTransferredMasses[pidi] + patchTransferredMassTotals[pidi]; + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace surfaceFilmModels +} // End namespace regionModels +} // End namespace Foam + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.H new file mode 100644 index 0000000000..f521672bb7 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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::regionModels::surfaceFilmModels::VoFPatchTransfer + +Description + Transfer mass between the film and the VoF in the continuous phase. + +SourceFiles + VoFPatchTransfer.C + +\*---------------------------------------------------------------------------*/ + +#ifndef VoFPatchTransfer_H +#define VoFPatchTransfer_H + +#include "transferModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace regionModels +{ +namespace surfaceFilmModels +{ + +/*---------------------------------------------------------------------------*\ + Class VoFPatchTransfer Declaration +\*---------------------------------------------------------------------------*/ + +class VoFPatchTransfer +: + public transferModel +{ + // Private member functions + + //- Disallow default bitwise copy construct + VoFPatchTransfer(const VoFPatchTransfer&); + + //- Disallow default bitwise assignment + void operator=(const VoFPatchTransfer&); + + +protected: + + //- Factor of the cell height above which the film is transferred + // to the VoF + scalar deltaFactorToVoF_; + + //- Factor of the cell height below which the VoF may be transferred + // to the film + scalar deltaFactorToFilm_; + + //- VoF limit above which all of the film is transferred to the VoF + scalar alphaToVoF_; + + //- VoF limit below which the VoF may be transferred to the film + scalar alphaToFilm_; + + //- Transfer rate coefficient + scalar transferRateCoeff_; + + //- List of patch IDs at which the film is removed + labelList patchIDs_; + + //- Transferred mass for each patch at which the film is removed + scalarField patchTransferredMasses_; + + +public: + + //- Runtime type information + TypeName("VoFPatchTransfer"); + + + // Constructors + + //- Construct from surface film model + VoFPatchTransfer(surfaceFilmRegionModel& film, const dictionary& dict); + + + //- Destructor + virtual ~VoFPatchTransfer(); + + + // Member Functions + + //- Correct + virtual void correct + ( + scalarField& availableMass, + scalarField& massToTransfer + ); + + //- Correct kinematic and thermodynamic transfers + virtual void correct + ( + scalarField& availableMass, + scalarField& massToTransfer, + scalarField& energyToTransfer + ); + + //- Accumulate the total mass injected for the patches into the + // scalarField provided + virtual void patchTransferredMassTotals + ( + scalarField& patchMasses + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace surfaceFilmModels +} // End namespace regionModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C new file mode 100644 index 0000000000..efcee3dca1 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C @@ -0,0 +1,215 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "VoFSolidificationMeltingSource.H" +#include "twoPhaseMixtureThermo.H" +#include "zeroGradientFvPatchFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +namespace Foam +{ + namespace fv + { + defineTypeNameAndDebug(VoFSolidificationMeltingSource, 0); + + addToRunTimeSelectionTable + ( + option, + VoFSolidificationMeltingSource, + dictionary + ); + } +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::fv::VoFSolidificationMeltingSource::update() +{ + if (curTimeIndex_ == mesh_.time().timeIndex()) + { + return; + } + + if (debug) + { + Info<< type() << ": " << name_ + << " - updating solid phase fraction" << endl; + } + + alphaSolid_.oldTime(); + + const twoPhaseMixtureThermo& thermo + ( + mesh_.lookupObject + ( + twoPhaseMixtureThermo::dictName + ) + ); + + const volScalarField& TVoF = thermo.thermo1().T(); + const volScalarField CpVoF(thermo.thermo1().Cp()); + const volScalarField& alphaVoF = thermo.alpha1(); + + forAll(cells_, i) + { + const label celli = cells_[i]; + + alphaSolid_[celli] = min + ( + relax_*alphaVoF[celli]*alphaSolidT_->value(TVoF[celli]) + + (1 - relax_)*alphaSolid_[celli], + alphaVoF[celli] + ); + } + + alphaSolid_.correctBoundaryConditions(); + + curTimeIndex_ = mesh_.time().timeIndex(); +} + + +Foam::word Foam::fv::VoFSolidificationMeltingSource::alphaSolidName() const +{ + const twoPhaseMixtureThermo& thermo + ( + mesh_.lookupObject + ( + twoPhaseMixtureThermo::dictName + ) + ); + + const volScalarField& alphaVoF = thermo.alpha1(); + + return IOobject::groupName(alphaVoF.name(), "solid"); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource +( + const word& sourceName, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh +) +: + cellSetOption(sourceName, modelType, dict, mesh), + alphaSolidT_(Function1::New("alphaSolidT", coeffs_)), + L_("L", dimEnergy/dimMass, coeffs_), + relax_(coeffs_.lookupOrDefault("relax", 0.9)), + Cu_(coeffs_.lookupOrDefault("Cu", 100000)), + q_(coeffs_.lookupOrDefault("q", 0.001)), + alphaSolid_ + ( + IOobject + ( + alphaSolidName(), + mesh.time().timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("alpha1", dimless, 0.0), + zeroGradientFvPatchScalarField::typeName + ), + curTimeIndex_(-1) +{ + fieldNames_.setSize(2); + fieldNames_[0] = "U"; + fieldNames_[1] = "T"; + applied_.setSize(fieldNames_.size(), false); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::fv::VoFSolidificationMeltingSource::addSup +( + fvMatrix& eqn, + const label fieldi +) +{ + apply(geometricOneField(), eqn); +} + + +void Foam::fv::VoFSolidificationMeltingSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi +) +{ + apply(rho, eqn); +} + + +void Foam::fv::VoFSolidificationMeltingSource::addSup +( + fvMatrix& eqn, + const label fieldi +) +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + update(); + + scalarField& Sp = eqn.diag(); + const scalarField& V = mesh_.V(); + + forAll(cells_, i) + { + const label celli = cells_[i]; + const scalar Vc = V[celli]; + const scalar alphaFluid = 1 - alphaSolid_[celli]; + + const scalar S = Cu_*sqr(1 - alphaFluid)/(pow3(alphaFluid) + q_); + + Sp[celli] -= Vc*S; + } +} + + +void Foam::fv::VoFSolidificationMeltingSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi +) +{ + // Momentum source uses a Boussinesq approximation - redirect + addSup(eqn, fieldi); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H new file mode 100644 index 0000000000..2749528e59 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H @@ -0,0 +1,215 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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::fv::VoFSolidificationMeltingSource + +Description + Solidification and melting model for VoF simulations. + + The presence of the solid phase in the flow field is incorporated into the + model as a momentum porosity contribution; the energy associated with the + phase change is added as an enthalpy contribution. The solid fraction as a + function of temperature \c alphaSolidT is specified as a Foam::Function1. + + The model writes the field \c alpha[01].solid which can be visualised to to + show the solid distribution. + +Usage + Example usage: + \verbatim + VoFSolidificationMeltingSource1 + { + type VoFSolidificationMeltingSource; + active yes; + + selectionMode cellZone; + cellZone solidZone; + + alphaSolidT table + ( + (330 1) + (335 0) + ); + + L 334000; + } + \endverbatim + + Where: + \table + Property | Description | Required | Default value + alphaSolidT | Solid fraction as function of temperature | yes | + L | Latent heat of fusion [J/kg] | yes | + relax | Relaxation coefficient [0-1] | no | 0.9 + Cu | Model coefficient | no | 100000 + q | Model coefficient | no | 0.001 + \endtable + +See also + Foam::fv::solidificationMeltingSource + Foam::Function1 + +SourceFiles + VoFSolidificationMeltingSource.C + VoFSolidificationMeltingSourceIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef VoFSolidificationMeltingSource_H +#define VoFSolidificationMeltingSource_H + +#include "fvMesh.H" +#include "volFields.H" +#include "cellSetOption.H" +#include "Function1.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class VoFSolidificationMeltingSource Declaration +\*---------------------------------------------------------------------------*/ + +class VoFSolidificationMeltingSource +: + public cellSetOption +{ + // Private data + + //- Solid fraction as a function of temperature + autoPtr> alphaSolidT_; + + //- Latent heat of fusion [J/kg] + dimensionedScalar L_; + + //- Phase fraction under-relaxation coefficient + scalar relax_; + + //- Mushy region momentum sink coefficient [1/s]; default = 10^5 + scalar Cu_; + + //- Coefficient used in porosity calc - default = 0.001 + scalar q_; + + //- Solid phase fraction + volScalarField alphaSolid_; + + //- Current time index (used for updating) + label curTimeIndex_; + + + // Private Member Functions + + //- Return the name of the solid phase fraction + word alphaSolidName() const; + + //- Update the model + void update(); + + //- Helper function to apply to the energy equation + template + void apply(const RhoFieldType& rho, fvMatrix& eqn); + + //- Disallow default bitwise copy construct + VoFSolidificationMeltingSource(const VoFSolidificationMeltingSource&); + + //- Disallow default bitwise assignment + void operator=(const VoFSolidificationMeltingSource&); + + +public: + + //- Runtime type information + TypeName("VoFSolidificationMeltingSource"); + + + // Constructors + + //- Construct from explicit source name and mesh + VoFSolidificationMeltingSource + ( + const word& sourceName, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + + // Member Functions + + // Add explicit and implicit contributions + + //- Add explicit contribution to enthalpy equation + virtual void addSup(fvMatrix& eqn, const label fieldi); + + //- Add implicit contribution to momentum equation + virtual void addSup(fvMatrix& eqn, const label fieldi); + + + // Add explicit and implicit contributions to compressible equation + + //- Add explicit contribution to compressible enthalpy equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi + ); + + //- Add implicit contribution to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldi + ); + + + // IO + + //- Read source dictionary + virtual bool read(const dictionary& dict); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "VoFSolidificationMeltingSourceTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C new file mode 100644 index 0000000000..43e5df3da0 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "VoFSolidificationMeltingSource.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::fv::VoFSolidificationMeltingSource::read(const dictionary& dict) +{ + if (cellSetOption::read(dict)) + { + alphaSolidT_ = Function1::New("alphaSolidT", coeffs_); + coeffs_.lookup("L") >> L_; + coeffs_.readIfPresent("relax", relax_); + coeffs_.readIfPresent("Cu", Cu_); + coeffs_.readIfPresent("q", q_); + + return true; + } + else + { + return false; + } + + return false; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C new file mode 100644 index 0000000000..a018696758 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C @@ -0,0 +1,66 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "fvcDdt.H" +#include "twoPhaseMixtureThermo.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::fv::VoFSolidificationMeltingSource::apply +( + const RhoFieldType& rho, + fvMatrix& eqn +) +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + update(); + + const twoPhaseMixtureThermo& thermo + ( + mesh_.lookupObject + ( + twoPhaseMixtureThermo::dictName + ) + ); + + const volScalarField CpVoF(thermo.thermo1().Cp()); + + if (eqn.psi().dimensions() == dimTemperature) + { + eqn += L_/CpVoF*(fvc::ddt(rho, alphaSolid_)); + } + else + { + eqn += L_*(fvc::ddt(rho, alphaSolid_)); + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C new file mode 100644 index 0000000000..bc30ed8f60 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2017 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 + compressibleInterFoam + +Description + 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 + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "CMULES.H" +#include "EulerDdtScheme.H" +#include "localEulerDdtScheme.H" +#include "CrankNicolsonDdtScheme.H" +#include "subCycle.H" +#include "rhoThermo.H" +#include "twoPhaseMixture.H" +#include "twoPhaseMixtureThermo.H" +#include "turbulentFluidThermoModel.H" +#include "SLGThermo.H" +#include "surfaceFilmModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "createFields.H" + #include "createAlphaFluxes.H" + #include "createSurfaceFilmModel.H" + #include "createFvOptions.H" + + volScalarField& p = mixture.p(); + volScalarField& T = mixture.T(); + const volScalarField& psi1 = mixture.thermo1().psi(); + const volScalarField& psi2 = mixture.thermo2().psi(); + + regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm(); + + turbulence->validate(); + + if (!LTS) + { + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + } + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + if (LTS) + { + #include "setRDeltaT.H" + } + else + { + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" + } + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + surfaceFilm.evolve(); + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "alphaControls.H" + #include "compressibleAlphaEqnSubCycle.H" + + volScalarField::Internal Srho(surfaceFilm.Srho()); + contErr -= posPart(Srho); + + #include "UEqn.H" + #include "TEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " + << runTime.elapsedCpuTime() + << " s\n\n" << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/createSurfaceFilmModel.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/createSurfaceFilmModel.H new file mode 100644 index 0000000000..1d452e7549 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/createSurfaceFilmModel.H @@ -0,0 +1,7 @@ +Info<< "\nConstructing surface film model" << endl; + +SLGThermo slgThermo(mesh, mixture.thermo1()); +autoPtr tsurfaceFilm +( + regionModels::surfaceFilmModel::New(mesh, g) +); diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H new file mode 100644 index 0000000000..27196f0eec --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H @@ -0,0 +1,129 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::flux(HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + MRF.makeRelative(phiHbyA); + + surfaceScalarField phig + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF); + + tmp p_rghEqnComp1; + tmp p_rghEqnComp2; + + if (pimple.transonic()) + { + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); + surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); + + p_rghEqnComp1 = + fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1) + + correction + ( + psi1*fvm::ddt(p_rgh) + + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + ); + deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr()); + p_rghEqnComp1.ref().relax(); + + p_rghEqnComp2 = + fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2) + + correction + ( + psi2*fvm::ddt(p_rgh) + + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) + ); + deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr()); + p_rghEqnComp2.ref().relax(); + } + else + { + #include "rhofs.H" + + p_rghEqnComp1 = + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)) + ) - surfaceFilm.Srho()/rho1; + + p_rghEqnComp2 = + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)) + ); + } + + // Cache p_rgh prior to solve for density update + volScalarField p_rgh_0(p_rgh); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqnIncomp + ( + fvc::div(phiHbyA) + - fvm::laplacian(rAUf, p_rgh) + ); + + solve + ( + p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp, + mesh.solver(p_rgh.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); + p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; + + dgdt = + ( + alpha1*(p_rghEqnComp2 & p_rgh) + - alpha2*(p_rghEqnComp1 & p_rgh) + ); + + phi = phiHbyA + p_rghEqnIncomp.flux(); + + U = HbyA + + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf); + U.correctBoundaryConditions(); + fvOptions.correct(U); + } + } + + // Update densities from change in p_rgh + mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0)); + mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0)); + + rho = alpha1*rho1 + alpha2*rho2; + + // Correct p_rgh for consistency with p and the updated densities + p_rgh = p - rho*gh; + p_rgh.correctBoundaryConditions(); + + K = 0.5*magSqr(U); +} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index dee1a578f2..e0e3c4459c 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C @@ -108,9 +108,7 @@ int main(int argc, char *argv[]) while (pimple.loop()) { #include "alphaControls.H" - #include "alphaEqnSubCycle.H" - - solve(fvm::ddt(rho) + fvc::div(rhoPhi)); + #include "compressibleAlphaEqnSubCycle.H" #include "UEqn.H" volScalarField divU(fvc::div(fvc::absolute(phi, U))); diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H index c1f3538169..8170ff2f99 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H @@ -87,10 +87,7 @@ surfaceScalarField rhoPhi fvc::interpolate(rho)*phi ); -volScalarField dgdt -( - pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001)) -); +volScalarField dgdt(alpha1*fvc::div(phi)); // Construct compressible turbulence model autoPtr turbulence diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 93412c9dba..54bbcc5e4d 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -53,13 +53,29 @@ } else { + #include "rhofs.H" + p_rghEqnComp1 = - fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh)) - + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1); + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)) + ); p_rghEqnComp2 = - fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh)) - + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2); + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)) + ); } // Cache p_rgh prior to solve for density update @@ -75,11 +91,7 @@ solve ( - ( - (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1() - + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2() - ) - + p_rghEqnIncomp, + p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp, mesh.solver(p_rgh.select(pimple.finalInnerIter())) ); @@ -90,8 +102,8 @@ dgdt = ( - pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 - - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 + alpha1*(p_rghEqnComp2 & p_rgh) + - alpha2*(p_rghEqnComp1 & p_rgh) ); phi = phiHbyA + p_rghEqnIncomp.flux(); diff --git a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C index 6274ad1c31..063db62589 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C @@ -63,8 +63,8 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo thermo1_ = rhoThermo::New(U.mesh(), phase1Name()); thermo2_ = rhoThermo::New(U.mesh(), phase2Name()); - thermo1_->validate(phase1Name(), "e"); - thermo2_->validate(phase2Name(), "e"); + // thermo1_->validate(phase1Name(), "e"); + // thermo2_->validate(phase2Name(), "e"); correct(); } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C index 9b064b14c6..5532f7a3d6 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -950,7 +950,7 @@ Foam::multiphaseMixtureThermo::nearInterface() const forAllConstIter(PtrDictionary, phases_, phase) { tnearInt.ref() = - max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase())); + max(tnearInt(), pos0(phase() - 0.01)*pos0(0.99 - phase())); } return tnearInt; diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H index 6526dacb58..f0d7e05ed2 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H @@ -104,8 +104,10 @@ ) { phase().dgdt() = - pos(phase()) + pos0(phase()) *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho(); + + phasei++; } phi = phiHbyA + p_rghEqnIncomp.flux(); diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 1dd92f95a1..615671b1a8 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -132,7 +132,7 @@ int main(int argc, char *argv[]) // if the mesh topology changed if (mesh.topoChanging()) { - talphaPhiCorr0.clear(); + talphaPhi1Corr0.clear(); } gh = (g & mesh.C()) - ghRef; diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H index 3f75765dd7..004c0d1ddb 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H @@ -183,8 +183,8 @@ solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1)); // Create the diffusion coefficients for alpha2<->alpha3 - volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2)); - volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3)); + volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2)); + volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3)); // Add the diffusive flux for alpha3->alpha2 alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1); diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C index b5495dd6ae..99b26c4a83 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C @@ -214,8 +214,8 @@ Foam::threePhaseInterfaceProperties::nearInterface() const { return max ( - pos(mixture_.alpha1() - 0.01)*pos(0.99 - mixture_.alpha1()), - pos(mixture_.alpha2() - 0.01)*pos(0.99 - mixture_.alpha2()) + pos0(mixture_.alpha1() - 0.01)*pos0(0.99 - mixture_.alpha1()), + pos0(mixture_.alpha2() - 0.01)*pos0(0.99 - mixture_.alpha2()) ); } diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C index 4a986ef175..c60310034e 100644 --- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C @@ -141,7 +141,7 @@ int main(int argc, char *argv[]) // if the mesh topology changed if (mesh.topoChanging()) { - talphaPhiCorr0.clear(); + talphaPhi1Corr0.clear(); } gh = (g & mesh.C()) - ghRef; diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C index 32d3fe4308..0f0cfac210 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C @@ -87,7 +87,7 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const return Pair> ( mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1) - *pos(p - pSat())/max(p - pSat(), 0.01*pSat()), + *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()), (-mvCoeff_)*limitedAlpha1*neg(p - pSat()) ); diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C index ade26e86a4..e9b1b93588 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C @@ -83,7 +83,7 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const return Pair> ( - mcCoeff_*(1.0 - limitedAlpha1)*pos(p - pSat()), + mcCoeff_*(1.0 - limitedAlpha1)*pos0(p - pSat()), (-mvCoeff_)*limitedAlpha1*neg(p - pSat()) ); } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C index 07fe4a2331..cc21a86a4f 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C @@ -136,7 +136,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const return Pair> ( - Cc_*(1.0 - limitedAlpha1)*pos(p - pSat())*apCoeff, + Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff, (-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff ); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index cc4c7f9656..cd9e68d7f8 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,13 +78,13 @@ Foam::tmp Foam::dragModels::GidaspowErgunWenYu::K volScalarField Cds ( neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) - + pos(Re - 1000)*0.44 + + pos0(Re - 1000)*0.44 ); // Wen and Yu (1966) return ( - pos(alpha2 - 0.8) + pos0(alpha2 - 0.8) *(0.75*Cds*phase2_.rho()*Ur*bp/d) + neg(alpha2 - 0.8) *( diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index 5e2048aad1..59d3358d4a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,7 +77,7 @@ Foam::tmp Foam::dragModels::GidaspowSchillerNaumann::K volScalarField Cds ( neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) - + pos(Re - 1000)*0.44 + + pos0(Re - 1000)*0.44 ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index d4d77a2bb3..a9983ed21a 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.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-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ Foam::tmp Foam::dragModels::SchillerNaumann::K volScalarField Cds ( neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) - + pos(Re - 1000)*0.44 + + pos0(Re - 1000)*0.44 ); return 0.75*Cds*phase2_.rho()*Ur/phase1_.d(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index 83419c54d1..975ad0c029 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -75,7 +75,7 @@ Foam::tmp Foam::dragModels::SyamlalOBrien::K volScalarField B ( neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28)) - + pos(alpha2 - 0.85)*(pow(alpha2, 2.65)) + + pos0(alpha2 - 0.85)*(pow(alpha2, 2.65)) ); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index fd644f5c5a..55c58b6c42 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,7 +77,7 @@ Foam::tmp Foam::dragModels::WenYu::K volScalarField Cds ( neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) - + pos(Re - 1000)*0.44 + + pos0(Re - 1000)*0.44 ); return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 425be38063..2f923c38a6 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -846,7 +846,8 @@ Foam::multiphaseSystem::nearInterface() const forAllConstIter(PtrDictionary, phases_, iter) { - tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter())); + tnearInt.ref() = + max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter())); } return tnearInt; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index 8e2d78cbec..c81f20db76 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -254,8 +254,8 @@ // dgdt = // ( - // pos(alpha2)*(pEqnComp2 & p)/rho2 - // - pos(alpha1)*(pEqnComp1 & p)/rho1 + // pos0(alpha2)*(pEqnComp2 & p)/rho2 + // - pos0(alpha1)*(pEqnComp1 & p)/rho1 // ); p_rgh.relax(); diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index 63a2776db4..825022dee9 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C @@ -548,7 +548,8 @@ Foam::multiphaseMixture::nearInterface() const forAllConstIter(PtrDictionary, phases_, iter) { - tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter())); + tnearInt.ref() = + max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter())); } return tnearInt; diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/files b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/files index 0642107d64..34919abb2b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/files +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/files @@ -17,6 +17,7 @@ saturationModels/Antoine/Antoine.C saturationModels/AntoineExtended/AntoineExtended.C saturationModels/ArdenBuck/ArdenBuck.C saturationModels/polynomial/polynomial.C +saturationModels/function1/function1.C saturationModels/constantSaturationConditions/constantSaturationConditions.C LIB = $(FOAM_LIBBIN)/libreactingEulerianInterfacialCompositionModels diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.C new file mode 100644 index 0000000000..f165c24d29 --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.C @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "function1.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace saturationModels +{ + defineTypeNameAndDebug(function1, 0); + addToRunTimeSelectionTable(saturationModel, function1, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::saturationModels::function1::function1(const dictionary& dict) +: + saturationModel(), + function_ + ( + Function1::New("function", dict) + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::saturationModels::function1::~function1() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::saturationModels::function1::pSat +( + const volScalarField& T +) const +{ + NotImplemented; + return volScalarField::null(); +} + + +Foam::tmp +Foam::saturationModels::function1::pSatPrime +( + const volScalarField& T +) const +{ + NotImplemented; + return volScalarField::null(); +} + + +Foam::tmp +Foam::saturationModels::function1::lnPSat +( + const volScalarField& T +) const +{ + NotImplemented; + return volScalarField::null(); +} + + +Foam::tmp +Foam::saturationModels::function1::Tsat +( + const volScalarField& p +) const +{ + tmp tTsat + ( + new volScalarField + ( + IOobject + ( + "Tsat", + p.mesh().time().timeName(), + p.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + p.mesh(), + dimensionedScalar("zero", dimTemperature, 0) + ) + ); + + volScalarField& Tsat = tTsat.ref(); + + forAll(Tsat, celli) + { + Tsat[celli] = function_->value(p[celli]); + } + + volScalarField::Boundary& TsatBf = Tsat.boundaryFieldRef(); + + forAll(Tsat.boundaryField(), patchi) + { + scalarField& Tsatp = TsatBf[patchi]; + const scalarField& pp = p.boundaryField()[patchi]; + + forAll(Tsatp, facei) + { + Tsatp[facei] = function_->value(pp[facei]); + + } + } + + return tTsat; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.H new file mode 100644 index 0000000000..bbac62d50b --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.H @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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::saturationModels::function1 + +Description + Saturation vapour temperature in terms of + the vapour pressure (in Pa). The saturation temperature in Kelvins is + specified as a Foam::Function1 type, to enable use of, e.g. constant, + polynomial, table values. + + Currently this class only provides \f$T_sat\f$, the inverse function to + return the vapour pressure for a given temperature are not implemented. + + Examples: + + \verbatim + type function1; + function polynomial + ( + (308.0422 0) + (0.0015096 1) + (-1.61589e-8 2) + (1.114106e-13 3) + (-4.52216e-19 4) + (1.05192e-24 5) + (-1.2953e-30 6) + (6.5365e-37 7) + ) + \endverbatim + + \verbatim + type function1; + + function csvFile; + functionCoeffs + { + nHeaderLine 1; + refColumn 0; + componentColumns (1); + separator ","; + mergeSeparators no; + file "filename.csv"; + outOfBounds clamp; + interpolationScheme linear; + }; + \endverbatim + +SourceFiles + function1.C + +\*---------------------------------------------------------------------------*/ + +#ifndef function1_saturationModel_H +#define function1_saturationModel_H + +#include "saturationModel.H" +#include "Function1.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace saturationModels +{ + +/*---------------------------------------------------------------------------*\ + Class function1 Declaration +\*---------------------------------------------------------------------------*/ + +class function1 +: + public saturationModel +{ + // Private data + + //- Saturation temperature as a function of pressure + autoPtr> function_; + + +public: + + //- Runtime type information + TypeName("function1"); + + // Constructors + + //- Construct from a dictionary + function1(const dictionary& dict); + + + //- Destructor + virtual ~function1(); + + + // Member Functions + + //- Saturation pressure + virtual tmp pSat(const volScalarField& T) const; + + //- Saturation pressure derivetive w.r.t. temperature + virtual tmp pSatPrime(const volScalarField& T) const; + + //- Natural log of the saturation pressure + virtual tmp lnPSat(const volScalarField& T) const; + + //- Saturation temperature + virtual tmp Tsat(const volScalarField& p) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace saturationModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C index 270419c5b4..073f70a5c7 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,9 +71,9 @@ Foam::aspectRatioModels::VakhrushevEfremov::E() const return neg(Ta - scalar(1))*scalar(1) - + pos(Ta - scalar(1))*neg(Ta - scalar(39.8)) + + pos0(Ta - scalar(1))*neg(Ta - scalar(39.8)) *pow3(0.81 + 0.206*tanh(1.6 - 2*log10(max(Ta, scalar(1))))) - + pos(Ta - scalar(39.8))*0.24; + + pos0(Ta - scalar(39.8))*0.24; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Beetstra/Beetstra.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Beetstra/Beetstra.C index 932bff3fb6..a76814e711 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Beetstra/Beetstra.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Beetstra/Beetstra.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -73,11 +73,12 @@ Foam::tmp Foam::dragModels::Beetstra::CdRe() const max(scalar(1) - pair_.dispersed(), pair_.continuous().residualAlpha()) ); - volScalarField Re(pair_.Re()); - volScalarField ReLim + volScalarField Res(alpha2*pair_.Re()); + + volScalarField ResLim ( "ReLim", - max(Re, residualRe_) + max(Res, residualRe_) ); volScalarField F0 @@ -89,9 +90,9 @@ Foam::tmp Foam::dragModels::Beetstra::CdRe() const volScalarField F1 ( "F1", - 0.413*Re/(24.0*sqr(alpha2))*(1.0/alpha2 - + 3.0*alpha1*alpha2 + 8.4*pow(ReLim, -0.343)) - /(1.0 + pow(10.0, 3*alpha1)*pow(ReLim, -(1.0 + 4.0*alpha1)/2.0)) + 0.413*Res/(24.0*sqr(alpha2))*(1.0/alpha2 + + 3.0*alpha1*alpha2 + 8.4*pow(ResLim, -0.343)) + /(1.0 + pow(10.0, 3*alpha1)*pow(ResLim, -(1.0 + 4.0*alpha1)/2.0)) ); return 24.0*alpha2*(F0 + F1); diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index 7649f217c8..6fec92d0b2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -84,7 +84,7 @@ Foam::tmp Foam::dragModels::GidaspowErgunWenYu::CdRe() const { return - pos(pair_.continuous() - 0.8)*WenYu_->CdRe() + pos0(pair_.continuous() - 0.8)*WenYu_->CdRe() + neg(pair_.continuous() - 0.8)*Ergun_->CdRe(); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index 9503264de9..792fc2a65e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ Foam::dragModels::GidaspowSchillerNaumann::CdRe() const volScalarField CdsRe ( neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))/alpha2 - + pos(Re - 1000)*0.44*max(Re, residualRe_) + + pos0(Re - 1000)*0.44*max(Re, residualRe_) ); return diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C index 8f6b115bbc..74adfb84a9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,7 +80,7 @@ Foam::dragModels::IshiiZuber::CdRe() const volScalarField ReM(Re*muc/muMix); volScalarField CdRe ( - pos(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687)) + pos0(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687)) + neg(1000 - ReM)*0.44*ReM ); @@ -92,7 +92,7 @@ Foam::dragModels::IshiiZuber::CdRe() const volScalarField CdReEllipse(Ealpha*0.6666*sqrt(Eo)*Re); return - pos(CdReEllipse - CdRe) + pos0(CdReEllipse - CdRe) *min(CdReEllipse, Re*sqr(1 - pair_.dispersed())*2.66667) + neg(CdReEllipse - CdRe)*CdRe; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Lain/Lain.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Lain/Lain.C index 2093120190..e172ecd7b2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Lain/Lain.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Lain/Lain.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,9 +66,9 @@ Foam::tmp Foam::dragModels::Lain::CdRe() const return neg(Re - 1.5)*16.0 - + pos(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22) - + pos(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL))) - + pos(Re - 1500.0)*2.61*Re; + + pos0(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22) + + pos0(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL))) + + pos0(Re - 1500.0)*2.61*Re; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index d0818e2c15..4ca93022f4 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,7 +67,7 @@ Foam::tmp Foam::dragModels::SchillerNaumann::CdRe() const return neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687)) - + pos(Re - 1000)*0.44*max(Re, residualRe_); + + pos0(Re - 1000)*0.44*max(Re, residualRe_); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index b722736b60..0c8c6f7b33 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,7 +71,7 @@ Foam::tmp Foam::dragModels::SyamlalOBrien::CdRe() const volScalarField B ( neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28)) - + pos(alpha2 - 0.85)*(pow(alpha2, 2.65)) + + pos0(alpha2 - 0.85)*(pow(alpha2, 2.65)) ); volScalarField Re(pair_.Re()); volScalarField Vr diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.C index 7fcda489a9..dcd0bba00e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,15 +50,6 @@ Foam::dragModels::Tenneti::Tenneti ) : dragModel(dict, pair, registerObject), - SchillerNaumann_ - ( - new SchillerNaumann - ( - dict, - pair, - false - ) - ), residualRe_("residualRe", dimless, dict.lookup("residualRe")) {} @@ -83,6 +74,14 @@ Foam::tmp Foam::dragModels::Tenneti::CdRe() const max(scalar(1) - pair_.dispersed(), pair_.continuous().residualAlpha()) ); + volScalarField Res(alpha2*pair_.Re()); + + volScalarField CdReIsolated + ( + neg(Res - 1000)*24.0*(1.0 + 0.15*pow(Res, 0.687)) + + pos0(Res - 1000)*0.44*max(Res, residualRe_) + ); + volScalarField F0 ( 5.81*alpha1/pow3(alpha2) + 0.48*pow(alpha1, 1.0/3.0)/pow4(alpha2) @@ -90,16 +89,14 @@ Foam::tmp Foam::dragModels::Tenneti::CdRe() const volScalarField F1 ( - pow(alpha1, 3)*max(pair_.Re(), residualRe_) - *(0.95 + 0.61*pow3(alpha1)/sqr(alpha2)) + pow3(alpha1)*Res*(0.95 + 0.61*pow3(alpha1)/sqr(alpha2)) ); // Tenneti et al. correlation includes the mean pressure drag. // This was removed here by multiplying F by alpha2 for consistency with // the formulation used in OpenFOAM return - SchillerNaumann_->CdRe()/(alpha2*max(pair_.Re(), residualRe_)) + - 24.0*sqr(alpha2)*(F0 + F1); + CdReIsolated + 24.0*sqr(alpha2)*(F0 + F1); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H index 7c9f64d57b..ae5c3c401f 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,9 +69,6 @@ class Tenneti { // Private data - //- SchillerNaumann drag model - autoPtr SchillerNaumann_; - //- Residual Reynolds Number const dimensionedScalar residualRe_; diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index c00679e5be..1d52dc6534 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,7 +72,7 @@ Foam::tmp Foam::dragModels::WenYu::CdRe() const volScalarField CdsRes ( neg(Res - 1000)*24.0*(1.0 + 0.15*pow(Res, 0.687)) - + pos(Res - 1000)*0.44*max(Res, residualRe_) + + pos0(Res - 1000)*0.44*max(Res, residualRe_) ); return diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C index d4a4dfd17e..ca541dbe41 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,8 +70,8 @@ Foam::tmp Foam::liftModels::TomiyamaLift::Cl() const return neg(EoH - scalar(4))*min(0.288*tanh(0.121*pair_.Re()), f) - + pos(EoH - scalar(4))*neg(EoH - scalar(10.7))*f - + pos(EoH - scalar(10.7))*(-0.288); + + pos0(EoH - scalar(4))*neg(EoH - scalar(10.7))*f + + pos0(EoH - scalar(10.7))*(-0.288); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C index 1f0fc837cb..d06a7929f2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -101,7 +101,11 @@ Foam::turbulentDispersionModels::Burns::D() const *sqr(pair_.dispersed().d()) ) *pair_.continuous().rho() - *(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_)); + *pair_.dispersed() + *( + 1.0/max(pair_.dispersed(), residualAlpha_) + + 1.0/max(pair_.continuous(), residualAlpha_) + ); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C index 6fc1940165..f549142810 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,9 +80,9 @@ Foam::tmp Foam::wallLubricationModels::Frank::Fi() const return zeroGradWalls ( ( - pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) - + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) - + pos(Eo - 33.0)*0.179 + pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) + + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) + + pos0(Eo - 33.0)*0.179 ) *max ( diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C index f6a1492900..47afcec18d 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,9 +78,9 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::Fi() const return zeroGradWalls ( ( - pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) - + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) - + pos(Eo - 33.0)*0.179 + pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) + + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) + + pos0(Eo - 33.0)*0.179 ) *0.5 *pair_.dispersed().d() diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index 4335bace8a..a6dd838d23 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -379,8 +379,8 @@ void Foam::ThermalPhaseChangePhaseSystem::correctThermo() ( min ( - (pos(iDmdt)*he2 + neg(iDmdt)*hef2) - - (neg(iDmdt)*he1 + pos(iDmdt)*hef1), + (pos0(iDmdt)*he2 + neg(iDmdt)*hef2) + - (neg(iDmdt)*he1 + pos0(iDmdt)*hef1), 0.3*mag(hef2 - hef1) ) ); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C index cabb9e2ebf..e7c8f0f8f9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -60,6 +60,13 @@ Foam::AnisothermalPhaseModel::~AnisothermalPhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +bool Foam::AnisothermalPhaseModel::compressible() const +{ + return !this->thermo().incompressible(); +} + + template void Foam::AnisothermalPhaseModel::correctKinematics() { @@ -135,7 +142,7 @@ Foam::AnisothermalPhaseModel::heEqn() he ) == - this->Qdot() + alpha*this->Qdot() ); // Add the appropriate pressure-work term @@ -156,13 +163,6 @@ Foam::AnisothermalPhaseModel::heEqn() } -template -bool Foam::AnisothermalPhaseModel::compressible() const -{ - return !this->thermo().incompressible(); -} - - template const Foam::volScalarField& Foam::AnisothermalPhaseModel::K() const diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H index 86edbcd1f8..95f7437f0a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,6 +85,9 @@ public: // Member Functions + //- Return true if the phase is compressible otherwise false + virtual bool compressible() const; + //- Correct the kinematics virtual void correctKinematics(); @@ -94,14 +97,8 @@ public: //- Return the enthalpy equation virtual tmp heEqn(); - - // Compressibility (variable density) - - //- Return true if the phase is compressible otherwise false - virtual bool compressible() const; - - //- Return the phase kinetic energy - virtual const volScalarField& K() const; + //- Return the phase kinetic energy + virtual const volScalarField& K() const; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C index 19b7b73ebe..95a9f0c38b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -49,6 +49,13 @@ Foam::IsothermalPhaseModel::~IsothermalPhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +bool Foam::IsothermalPhaseModel::compressible() const +{ + return !this->thermo().incompressible(); +} + + template void Foam::IsothermalPhaseModel::correctThermo() {} diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H index 85a31f8807..f8bb151fe3 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,6 +71,9 @@ public: // Member Functions + //- Return true if the phase is compressible otherwise false + virtual bool compressible() const; + //- Correct the thermodynamics virtual void correctThermo(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C index ba7af6a523..279b3bb64a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -172,7 +172,7 @@ Foam::MultiComponentPhaseModel::YiEqn Yi ) == - this->R(Yi) + alpha*this->R(Yi) + fvc::ddt(residualAlpha_*rho, Yi) - fvm::ddt(residualAlpha_*rho, Yi) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H index e4689efee6..dc53815cf4 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H @@ -31,6 +31,7 @@ for (int Ecorr=0; Ecorrrelax(); fvOptions.constrain(EEqn.ref()); EEqn->solve(); + fvOptions.correct(phase.thermo().he()); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 1f8f83599b..d9aeeded84 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -598,7 +598,7 @@ Foam::multiphaseSystem::nearInterface() const tnearInt.ref() = max ( tnearInt(), - pos(phases()[phasei] - 0.01)*pos(0.99 - phases()[phasei]) + pos0(phases()[phasei] - 0.01)*pos0(0.99 - phases()[phasei]) ); } @@ -699,6 +699,9 @@ void Foam::multiphaseSystem::solve() { phaseModel& phase = phases()[phasei]; phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi(); + + // Ensure the phase-fractions are bounded + phase.maxMin(0, 1); } calcAlphas(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index 3445c81fb9..a164b73f26 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -179,7 +179,7 @@ while (pimple.correct()) ( fvc::interpolate(fvc::average(alphafs[phasei])) ); - surfaceScalarField phiCorrCoeff(pos(alphafBar - 0.99)); + surfaceScalarField phiCorrCoeff(pos0(alphafBar - 0.99)); surfaceScalarField::Boundary& phiCorrCoeffBf = phiCorrCoeff.boundaryFieldRef(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H index d9f1e85e6f..246bea1fb4 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H @@ -25,6 +25,7 @@ for (int Ecorr=0; Ecorrrelax(); fvOptions.constrain(E1Eqn.ref()); E1Eqn->solve(); + fvOptions.correct(phase1.thermo().he()); } } @@ -45,6 +46,7 @@ for (int Ecorr=0; Ecorrrelax(); fvOptions.constrain(E2Eqn.ref()); E2Eqn->solve(); + fvOptions.correct(phase2.thermo().he()); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H index 3365febdf6..c30d4f1465 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -143,8 +143,8 @@ while (pimple.correct()) // ddtPhiCorr filter -- only apply in pure(ish) phases surfaceScalarField alphaf1Bar(fvc::interpolate(fvc::average(alphaf1))); - surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99)); - surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar)); + surfaceScalarField phiCorrCoeff1(pos0(alphaf1Bar - 0.99)); + surfaceScalarField phiCorrCoeff2(pos0(0.01 - alphaf1Bar)); { surfaceScalarField::Boundary& phiCorrCoeff1Bf = diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.C index 1963b11425..04b8c2917e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.C @@ -73,7 +73,7 @@ Lavieville::fLiquid ) const { return - pos(alphaLiquid - alphaCrit_) + pos0(alphaLiquid - alphaCrit_) *( 1 - 0.5*exp(-20*(alphaLiquid - alphaCrit_)) ) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C index 702464b813..1d157b67ab 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ cosine::fLiquid ) const { return - pos(alphaLiquid1_ - alphaLiquid) + pos0(alphaLiquid1_ - alphaLiquid) *( neg(alphaLiquid0_ - alphaLiquid) *( diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C index 2d8a5cce8e..0fa9d5faf6 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -267,7 +267,7 @@ void Foam::JohnsonJacksonParticleThetaFvPatchScalarField::updateCoeffs() this->refValue() = 0.0; this->refGrad() = - pos(alpha - SMALL) + pos0(alpha - SMALL) *constant::mathematical::pi *specularityCoefficient_.value() *alpha diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 62bbe6af40..1ab0d22780 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -403,8 +403,7 @@ void Foam::twoPhaseSystem::solve() << endl; // Ensure the phase-fractions are bounded - alpha1.max(0); - alpha1.min(1); + alpha1.maxMin(0, 1); // Update the phase-fraction of the other phase alpha2 = scalar(1) - alpha1; diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaCourantNo.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaCourantNo.H index 4e13f6ac56..9c8e27d975 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaCourantNo.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaCourantNo.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,7 +41,7 @@ if (mesh.nInternalFaces()) { scalarField sumPhi ( - pos(alpha1 - 0.01)*pos(0.99 - alpha1) + pos0(alpha1 - 0.01)*pos0(0.99 - alpha1) *fvc::surfaceSum(mag(phi))().primitiveField() ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H index 308e7589f9..e6dec13b3f 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H @@ -75,9 +75,11 @@ fvOptions.constrain(E1Eqn); E1Eqn.solve(); + fvOptions.correct(he1); fvOptions.constrain(E2Eqn); E2Eqn.solve(); + fvOptions.correct(he2); thermo1.correct(); Info<< "min " << thermo1.T().name() diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C index 270419c5b4..073f70a5c7 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,9 +71,9 @@ Foam::aspectRatioModels::VakhrushevEfremov::E() const return neg(Ta - scalar(1))*scalar(1) - + pos(Ta - scalar(1))*neg(Ta - scalar(39.8)) + + pos0(Ta - scalar(1))*neg(Ta - scalar(39.8)) *pow3(0.81 + 0.206*tanh(1.6 - 2*log10(max(Ta, scalar(1))))) - + pos(Ta - scalar(39.8))*0.24; + + pos0(Ta - scalar(39.8))*0.24; } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index ed8a871416..3ad9969954 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -84,7 +84,7 @@ Foam::tmp Foam::dragModels::GidaspowErgunWenYu::CdRe() const { return - pos(pair_.continuous() - 0.8)*WenYu_->CdRe() + pos0(pair_.continuous() - 0.8)*WenYu_->CdRe() + neg(pair_.continuous() - 0.8)*Ergun_->CdRe(); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index 9503264de9..792fc2a65e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ Foam::dragModels::GidaspowSchillerNaumann::CdRe() const volScalarField CdsRe ( neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))/alpha2 - + pos(Re - 1000)*0.44*max(Re, residualRe_) + + pos0(Re - 1000)*0.44*max(Re, residualRe_) ); return diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C index 18971fbdc9..74adfb84a9 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,7 +80,7 @@ Foam::dragModels::IshiiZuber::CdRe() const volScalarField ReM(Re*muc/muMix); volScalarField CdRe ( - pos(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687)) + pos0(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687)) + neg(1000 - ReM)*0.44*ReM ); @@ -92,7 +92,7 @@ Foam::dragModels::IshiiZuber::CdRe() const volScalarField CdReEllipse(Ealpha*0.6666*sqrt(Eo)*Re); return - pos(CdReEllipse - CdRe) + pos0(CdReEllipse - CdRe) *min(CdReEllipse, Re*sqr(1 - pair_.dispersed())*2.66667) + neg(CdReEllipse - CdRe)*CdRe; } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Lain/Lain.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Lain/Lain.C index 74ce22f7b2..e172ecd7b2 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Lain/Lain.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Lain/Lain.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,9 +66,9 @@ Foam::tmp Foam::dragModels::Lain::CdRe() const return neg(Re - 1.5)*16.0 - + pos(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22) - + pos(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL))) - + pos(Re - 1500.0)*2.61*Re; + + pos0(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22) + + pos0(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL))) + + pos0(Re - 1500.0)*2.61*Re; } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index d0818e2c15..4ca93022f4 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,7 +67,7 @@ Foam::tmp Foam::dragModels::SchillerNaumann::CdRe() const return neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687)) - + pos(Re - 1000)*0.44*max(Re, residualRe_); + + pos0(Re - 1000)*0.44*max(Re, residualRe_); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index b722736b60..0c8c6f7b33 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,7 +71,7 @@ Foam::tmp Foam::dragModels::SyamlalOBrien::CdRe() const volScalarField B ( neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28)) - + pos(alpha2 - 0.85)*(pow(alpha2, 2.65)) + + pos0(alpha2 - 0.85)*(pow(alpha2, 2.65)) ); volScalarField Re(pair_.Re()); volScalarField Vr diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index c00679e5be..1d52dc6534 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,7 +72,7 @@ Foam::tmp Foam::dragModels::WenYu::CdRe() const volScalarField CdsRes ( neg(Res - 1000)*24.0*(1.0 + 0.15*pow(Res, 0.687)) - + pos(Res - 1000)*0.44*max(Res, residualRe_) + + pos0(Res - 1000)*0.44*max(Res, residualRe_) ); return diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C index d4a4dfd17e..ca541dbe41 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,8 +70,8 @@ Foam::tmp Foam::liftModels::TomiyamaLift::Cl() const return neg(EoH - scalar(4))*min(0.288*tanh(0.121*pair_.Re()), f) - + pos(EoH - scalar(4))*neg(EoH - scalar(10.7))*f - + pos(EoH - scalar(10.7))*(-0.288); + + pos0(EoH - scalar(4))*neg(EoH - scalar(10.7))*f + + pos0(EoH - scalar(10.7))*(-0.288); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C index 093013cca3..33db7656cb 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -101,7 +101,11 @@ Foam::turbulentDispersionModels::Burns::D() const *sqr(pair_.dispersed().d()) ) *pair_.continuous().rho() - *(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_)); + *pair_.dispersed() + *( + 1.0/max(pair_.dispersed(), residualAlpha_) + + 1.0/max(pair_.continuous(), residualAlpha_) + ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C index 0b9065a1a9..024e5ef2df 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -79,9 +79,9 @@ Foam::tmp Foam::wallLubricationModels::Frank::Fi() const return ( - pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) - + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) - + pos(Eo - 33.0)*0.179 + pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) + + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) + + pos0(Eo - 33.0)*0.179 ) *max ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C index 9209c5108c..3abdfa3d4e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,9 +77,9 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::Fi() const return ( - pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) - + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) - + pos(Eo - 33.0)*0.179 + pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179) + + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187) + + pos0(Eo - 33.0)*0.179 ) *0.5 *pair_.dispersed().d() diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H index 2010e65ebc..eabc3961df 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H @@ -145,8 +145,8 @@ while (pimple.correct()) // ddtPhiCorr filter -- only apply in pure(ish) phases surfaceScalarField alphaf1Bar(fvc::interpolate(fvc::average(alphaf1))); - surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99)); - surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar)); + surfaceScalarField phiCorrCoeff1(pos0(alphaf1Bar - 0.99)); + surfaceScalarField phiCorrCoeff2(pos0(0.01 - alphaf1Bar)); { surfaceScalarField::Boundary& phiCorrCoeff1Bf = diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C index 2d8a5cce8e..0fa9d5faf6 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2014-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -267,7 +267,7 @@ void Foam::JohnsonJacksonParticleThetaFvPatchScalarField::updateCoeffs() this->refValue() = 0.0; this->refGrad() = - pos(alpha - SMALL) + pos0(alpha - SMALL) *constant::mathematical::pi *specularityCoefficient_.value() *alpha diff --git a/applications/test/CompactIOList/Test-CompactIOList.C b/applications/test/CompactIOList/Test-CompactIOList.C index 12d64c6bde..2024e663e0 100644 --- a/applications/test/CompactIOList/Test-CompactIOList.C +++ b/applications/test/CompactIOList/Test-CompactIOList.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -87,7 +87,8 @@ int main(int argc, char *argv[]) ( format, IOstream::currentVersion, - IOstream::UNCOMPRESSED + IOstream::UNCOMPRESSED, + true ); Info<< "Written old format faceList in = " @@ -149,7 +150,8 @@ int main(int argc, char *argv[]) ( format, IOstream::currentVersion, - IOstream::UNCOMPRESSED + IOstream::UNCOMPRESSED, + true ); Info<< "Written new format faceList in = " diff --git a/applications/test/IOField/Make/files b/applications/test/IOField/Make/files new file mode 100644 index 0000000000..74e614deae --- /dev/null +++ b/applications/test/IOField/Make/files @@ -0,0 +1,3 @@ +Test-IOField.C + +EXE = $(FOAM_USER_APPBIN)/Test-IOField diff --git a/applications/test/IOField/Make/options b/applications/test/IOField/Make/options new file mode 100644 index 0000000000..6a9e9810b3 --- /dev/null +++ b/applications/test/IOField/Make/options @@ -0,0 +1,2 @@ +/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */ +/* EXE_LIBS = -lfiniteVolume */ diff --git a/applications/test/IOField/Test-IOField.C b/applications/test/IOField/Test-IOField.C new file mode 100644 index 0000000000..d49acaf8f5 --- /dev/null +++ b/applications/test/IOField/Test-IOField.C @@ -0,0 +1,189 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 + Test-IOField + +Description + Test the processor-local reading of IOField (used in the lagrangian libs) + +\*---------------------------------------------------------------------------*/ + +#include "IOField.H" +#include "argList.H" +#include "polyMesh.H" +#include "Time.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void write(const IOobject& io, const label sz) +{ + IOField