diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C index 0ddfa64588..000948514e 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C +++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C @@ -30,7 +30,7 @@ Description Sub-models include: - turbulence modelling, i.e. laminar, RAS or LES - - run-time selectable fvOptions, e.g. MRF, explicit porosity + - run-time selectable finitie volume options, e.g. MRF, explicit porosity \*---------------------------------------------------------------------------*/ diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Allwmake b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake index ecd3f26015..f615df216c 100755 --- a/applications/solvers/lagrangian/reactingParcelFoam/Allwmake +++ b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake @@ -3,6 +3,7 @@ cd ${0%/*} || exit 1 # run from this directory set -x wmake +wmake simpleReactingParcelFoam wmake LTSReactingParcelFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C index 991de15ad6..34aecb65b4 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C @@ -28,7 +28,7 @@ Description Local time stepping (LTS) solver for steady, compressible, laminar or turbulent reacting and non-reacting flow with multiphase Lagrangian parcels and porous media, including run-time selectable finitite volume - options, e.g. fvOptions, constraints + options, e.g. sources, constraints Note: ddtPhiCorr not used here when porous zones are active - not well defined for porous calculations diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C index 424c23b8b9..c4503cc101 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C @@ -22,12 +22,12 @@ License along with OpenFOAM. If not, see . Application - porousExplicitSourceReactingParcelFoam + reactingParcelFoam Description Transient PIMPLE solver for compressible, laminar or turbulent flow with reacting multiphase Lagrangian parcels, including run-time selectable - finite volume options, e.g. fvOptions, constraints + finite volume options, e.g. sources, constraints \*---------------------------------------------------------------------------*/ diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H new file mode 100644 index 0000000000..d4d686e355 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H @@ -0,0 +1,32 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + mvConvection->fvmDiv(phi, he) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho)) + : fvc::div(phi, volScalarField("K", 0.5*magSqr(U))) + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + parcels.Sh(he) + + radiation->Sh(thermo) + + combustion->Sh() + + fvOptions(rho, he) + ); + + EEqn.relax(); + + fvOptions.constrain(EEqn); + + EEqn.solve(); + + fvOptions.correct(he); + thermo.correct(); + radiation->correct(); + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; +} diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files new file mode 100644 index 0000000000..4a202fcd4d --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files @@ -0,0 +1,3 @@ +simpleReactingParcelFoam.C + +EXE = $(FOAM_APPBIN)/simpleReactingParcelFoam diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options new file mode 100644 index 0000000000..705cfc9433 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options @@ -0,0 +1,53 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I${LIB_SRC}/meshTools/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -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)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/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/radiationModels/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/combustionModels/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(FOAM_SOLVERS)/combustion/reactingFoam + + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -llagrangian \ + -llagrangianIntermediate \ + -lspecie \ + -lfluidThermophysicalModels \ + -lliquidProperties \ + -lliquidMixtureProperties \ + -lsolidProperties \ + -lsolidMixtureProperties \ + -lthermophysicalFunctions \ + -lreactionThermophysicalModels \ + -lSLGThermo \ + -lchemistryModel \ + -lradiationModels \ + -lODE \ + -lregionModels \ + -lsurfaceFilmModels \ + -lcombustionModels \ + -lfvOptions \ + -lsampling diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H new file mode 100644 index 0000000000..eabbb40454 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H @@ -0,0 +1,17 @@ + tmp UEqn + ( + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) + == + rho.dimensionedInternalField()*g + + parcels.SU(U) + + fvOptions(rho, U) + ); + + UEqn().relax(); + + fvOptions.constrain(UEqn()); + + solve(UEqn() == -fvc::grad(p)); + + fvOptions.correct(U); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H new file mode 100644 index 0000000000..cd0a45f0f0 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H @@ -0,0 +1,53 @@ +tmp > mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); + +{ + combustion->correct(); + dQ = combustion->dQ(); + label inertIndex = -1; + volScalarField Yt(0.0*Y[0]); + + forAll(Y, i) + { + if (Y[i].name() != inertSpecie) + { + volScalarField& Yi = Y[i]; + + fvScalarMatrix YEqn + ( + mvConvection->fvmDiv(phi, Yi) + - fvm::laplacian(turbulence->muEff(), Yi) + == + parcels.SYi(i, Yi) + + combustion->R(Yi) + + fvOptions(rho, Yi) + ); + + YEqn.relax(); + + fvOptions.constrain(YEqn); + + YEqn.solve(mesh.solver("Yi")); + + fvOptions.correct(Yi); + + Yi.max(0.0); + Yt += Yi; + } + else + { + inertIndex = i; + } + } + + Y[inertIndex] = scalar(1) - Yt; + Y[inertIndex].max(0.0); +} diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H new file mode 100644 index 0000000000..954b74e069 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H @@ -0,0 +1,9 @@ +Info<< "\nConstructing reacting cloud" << endl; +basicReactingMultiphaseCloud parcels +( + "reactingCloud1", + rho, + U, + g, + slgThermo +); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H new file mode 100644 index 0000000000..d6df24cb48 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H @@ -0,0 +1,98 @@ + Info<< "Creating combustion model\n" << endl; + + autoPtr combustion + ( + combustionModels::rhoCombustionModel::New(mesh) + ); + + rhoReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); + + SLGThermo slgThermo(mesh, thermo); + + basicMultiComponentMixture& composition = thermo.composition(); + PtrList& Y = composition.Y(); + + const word inertSpecie(thermo.lookup("inertSpecie")); + + if (!composition.contains(inertSpecie)) + { + FatalErrorIn(args.executable()) + << "Specified inert specie '" << inertSpecie << "' not found in " + << "species list. Available species:" << composition.species() + << exit(FatalError); + } + + volScalarField& p = thermo.p(); + const volScalarField& T = thermo.T(); + const volScalarField& psi = thermo.psi(); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + Info<< "\nReading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "compressibleCreatePhi.H" + + dimensionedScalar rhoMax(simple.dict().lookup("rhoMax")); + dimensionedScalar rhoMin(simple.dict().lookup("rhoMin")); + + 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()); + + Info<< "Creating multi-variate interpolation scheme\n" << endl; + multivariateSurfaceInterpolationScheme::fieldTable fields; + + forAll(Y, i) + { + fields.add(Y[i]); + } + fields.add(thermo.he()); + + volScalarField dQ + ( + IOobject + ( + "dQ", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) + ); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H new file mode 100644 index 0000000000..0d1aa2e238 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H @@ -0,0 +1,59 @@ +{ + rho = thermo.rho(); + + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + thermo.rho() -= psi*p; + + volScalarField rAU(1.0/UEqn().A()); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn().H(); + + UEqn.clear(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) + == + parcels.Srho() + + fvOptions(psi, p, rho.name()) + ); + + fvOptions.constrain(pEqn); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } + + p.relax(); + + // Second part of thermodynamic density update + thermo.rho() += psi*p; + + #include "compressibleContinuityErrs.H" + + U = HbyA - rAU*fvc::grad(p); + U.correctBoundaryConditions(); + fvOptions.correct(U); + + 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/simpleReactingParcelFoam/simpleReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C new file mode 100644 index 0000000000..6620d2af52 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 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 + simpleReactingParcelFoam + +Description + Steady state SIMPLE solver for compressible, laminar or turbulent flow with + reacting multiphase Lagrangian parcels, including run-time selectable + finite volume options, e.g. sources, constraints + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "turbulenceModel.H" +#include "basicReactingMultiphaseCloud.H" +#include "rhoCombustionModel.H" +#include "radiationModel.H" +#include "IOporosityModelList.H" +#include "fvIOoptionList.H" +#include "SLGThermo.H" +#include "simpleControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + + simpleControl simple(mesh); + + #include "createFields.H" + #include "createRadiationModel.H" + #include "createClouds.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + parcels.evolve(); + + // --- Pressure-velocity SIMPLE corrector loop + { + #include "UEqn.H" + #include "YEqn.H" + #include "EEqn.H" + #include "pEqn.H" + } + + turbulence->correct(); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H new file mode 100644 index 0000000000..c894b36b41 --- /dev/null +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H @@ -0,0 +1,49 @@ +{ + volScalarField k1 + ( + "k1", + alpha1*(thermo1.alpha()/rho1 + sqr(Ct)*nut2*thermo1.CpByCpv()/Prt) + ); + + volScalarField k2 + ( + "k2", + alpha2*(thermo2.alpha()/rho2 + nut2*thermo2.CpByCpv()/Prt) + ); + + volScalarField& he1 = thermo1.he(); + volScalarField& he2 = thermo2.he(); + + Info<< max(he1) << min(he1) << endl; + + fvScalarMatrix he1Eqn + ( + fvm::ddt(alpha1, he1) + + fvm::div(alphaPhi1, he1) + - fvm::laplacian(k1, he1) + == + heatTransferCoeff*(he1/thermo1.Cp())/rho1 + - fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1) + + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1))) + ); + + fvScalarMatrix he2Eqn + ( + fvm::ddt(alpha2, he2) + + fvm::div(alphaPhi2, he2) + - fvm::laplacian(k2, he2) + == + heatTransferCoeff*(he2/thermo2.Cp())/rho2 + - fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2) + + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2))) + ); + + he1Eqn.relax(); + he1Eqn.solve(); + + he2Eqn.relax(); + he2Eqn.solve(); + + thermo1.correct(); + thermo2.correct(); +} diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options index 16ce2e3674..d90ba7bbbd 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -IturbulenceModel \ @@ -8,6 +9,8 @@ EXE_INC = \ -Iaveraging EXE_LIBS = \ + -lfluidThermophysicalModels \ + -lspecie \ -lincompressibleTransportModels \ -lcompressiblePhaseModel \ -lcompressibleEulerianInterfacialModels \ diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H deleted file mode 100644 index 8f38ca9d91..0000000000 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H +++ /dev/null @@ -1,36 +0,0 @@ -{ - volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt)); - volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt)); - - fvScalarMatrix T1Eqn - ( - fvm::ddt(alpha1, T1) - + fvm::div(alphaPhi1, T1) - - fvm::laplacian(kByCp1, T1) - == - heatTransferCoeff*T2/Cp1/rho1 - - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1) - + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1 - ); - - fvScalarMatrix T2Eqn - ( - fvm::ddt(alpha2, T2) - + fvm::div(alphaPhi2, T2) - - fvm::laplacian(kByCp2, T2) - == - heatTransferCoeff*T1/Cp2/rho2 - - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2) - + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2 - ); - - T1Eqn.relax(); - T1Eqn.solve(); - - T2Eqn.relax(); - T2Eqn.solve(); - - // Update compressibilities - psi1 = 1.0/(R1*T1); - psi2 = 1.0/(R2*T2); -} diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H index 606b4573e4..a243db6300 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H @@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); } else // If not using kinetic theory is using Ct model { - nuEff1 = sqr(Ct)*nut2 + nu1; + nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1; } volTensorField Rc1 ( - "Rc1", + "Rc", (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T ); @@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); volTensorField gradU2T(fvc::grad(U2)().T()); volTensorField Rc2 ( - "Rc2", + "Rc", (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H index a818ee2e9e..86d9203dce 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H @@ -1,9 +1,9 @@ -surfaceScalarField alphaPhi1("alphaPhi1", phi1); -surfaceScalarField alphaPhi2("alphaPhi2", phi2); +surfaceScalarField alphaPhi1("alphaPhi", phi1); +surfaceScalarField alphaPhi2("alphaPhi", phi2); { - word scheme("div(phi,alpha1)"); - word schemer("div(phir,alpha1)"); + word scheme("div(phi,alpha)"); + word schemer("div(phir,alpha)"); surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); @@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2); surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf()); phir += phipp; - phic += fvc::interpolate(alpha1)*phipp; + phic += alpha1f*phipp; } for (int acorr=0; acorr 0.0) { + surfaceScalarField alpha1f(fvc::interpolate(alpha1)); + ppMagf = fvc::interpolate((1.0/rho1)*rAU1) - *fvc::interpolate - ( - g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) - ); + *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax); + + // ppMagf = + // fvc::interpolate((1.0/rho1)*rAU1) + // *fvc::interpolate + // ( + // g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) + // ); alpha1Eqn -= fvm::laplacian ( - (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, + alpha1f*ppMagf, alpha1, "laplacian(alphaPpMag,alpha1)" ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C index 855a862b6b..65372612b6 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -31,6 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "rhoThermo.H" #include "nearWallDist.H" #include "wallFvPatch.H" #include "fixedValueFvsPatchFields.H" @@ -86,7 +87,7 @@ int main(int argc, char *argv[]) #include "alphaEqn.H" #include "kEpsilon.H" #include "interfacialCoeffs.H" - #include "TEqns.H" + #include "EEqns.H" #include "UEqns.H" // --- Pressure corrector loop diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H index e0cd5328c9..6095ad6845 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H @@ -12,55 +12,43 @@ ) ); + word phase1Name + ( + transportProperties.found("phases") + ? wordList(transportProperties.lookup("phases"))[0] + : "phase1" + ); + + word phase2Name + ( + transportProperties.found("phases") + ? wordList(transportProperties.lookup("phases"))[1] + : "phase2" + ); + autoPtr phase1 = phaseModel::New ( mesh, transportProperties, - "1" + phase1Name ); autoPtr phase2 = phaseModel::New ( mesh, transportProperties, - "2" + phase2Name ); + volScalarField& alpha1 = phase1(); + volScalarField& alpha2 = phase2(); + alpha2 = scalar(1) - alpha1; + volVectorField& U1 = phase1->U(); surfaceScalarField& phi1 = phase1->phi(); - const dimensionedScalar& nu1 = phase1->nu(); - const dimensionedScalar& k1 = phase1->kappa(); - const dimensionedScalar& Cp1 = phase1->Cp(); - dimensionedScalar rho10 - ( - "rho0", - dimDensity, - transportProperties.subDict("phase1").lookup("rho0") - ); - dimensionedScalar R1 - ( - "R", - dimensionSet(0, 2, -2, -1, 0), - transportProperties.subDict("phase1").lookup("R") - ); volVectorField& U2 = phase2->U(); surfaceScalarField& phi2 = phase2->phi(); - const dimensionedScalar& nu2 = phase2->nu(); - const dimensionedScalar& k2 = phase2->kappa(); - const dimensionedScalar& Cp2 = phase2->Cp(); - dimensionedScalar rho20 - ( - "rho0", - dimDensity, - transportProperties.subDict("phase2").lookup("rho0") - ); - dimensionedScalar R2 - ( - "R", - dimensionSet(0, 2, -2, -1, 0), - transportProperties.subDict("phase2").lookup("R") - ); dimensionedScalar pMin ( @@ -69,102 +57,16 @@ transportProperties.lookup("pMin") ); + rhoThermo& thermo1 = phase1->thermo(); + rhoThermo& thermo2 = phase2->thermo(); - Info<< "Reading field T1" << endl; - volScalarField T1 - ( - IOobject - ( - "T1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); + volScalarField& p = thermo1.p(); - Info<< "Reading field T2" << endl; - volScalarField T2 - ( - IOobject - ( - "T2", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); + volScalarField rho1("rho" + phase1Name, thermo1.rho()); + const volScalarField& psi1 = thermo1.psi(); - volScalarField psi1 - ( - IOobject - ( - "psi1", - runTime.timeName(), - mesh - ), - 1.0/(R1*T1) - ); - - volScalarField psi2 - ( - IOobject - ( - "psi2", - runTime.timeName(), - mesh - ), - 1.0/(R2*T2) - ); - psi2.oldTime(); - - Info<< "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField rho1("rho1", rho10 + psi1*p); - volScalarField rho2("rho2", rho20 + psi2*p); - - - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField alpha2 - ( - IOobject - ( - "alpha2", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - scalar(1) - alpha1 - ); + volScalarField rho2("rho" + phase2Name, thermo2.rho()); + const volScalarField& psi2 = thermo2.psi(); volVectorField U ( @@ -179,27 +81,6 @@ alpha1*U1 + alpha2*U2 ); - dimensionedScalar Cvm - ( - "Cvm", - dimless, - transportProperties.lookup("Cvm") - ); - - dimensionedScalar Cl - ( - "Cl", - dimless, - transportProperties.lookup("Cl") - ); - - dimensionedScalar Ct - ( - "Ct", - dimless, - transportProperties.lookup("Ct") - ); - surfaceScalarField phi ( IOobject @@ -226,8 +107,6 @@ alpha1*rho1 + alpha2*rho2 ); - #include "createRASTurbulence.H" - Info<< "Calculating field DDtU1 and DDtU2\n" << endl; volVectorField DDtU1 @@ -248,6 +127,29 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + dimensionedScalar Cvm + ( + "Cvm", + dimless, + transportProperties.lookup("Cvm") + ); + + dimensionedScalar Cl + ( + "Cl", + dimless, + transportProperties.lookup("Cl") + ); + + dimensionedScalar Ct + ( + "Ct", + dimless, + transportProperties.lookup("Ct") + ); + + #include "createRASTurbulence.H" + IOdictionary interfacialProperties ( IOobject @@ -297,8 +199,8 @@ if ( !( - dispersedPhase == "1" - || dispersedPhase == "2" + dispersedPhase == phase1Name + || dispersedPhase == phase2Name || dispersedPhase == "both" ) ) @@ -337,7 +239,7 @@ ( IOobject ( - "rAU1", + "rAU" + phase1Name, runTime.timeName(), mesh, IOobject::NO_READ, @@ -387,5 +289,5 @@ ); Info<< "Creating field kinetic energy K\n" << endl; - volScalarField K1("K1", 0.5*magSqr(U1)); - volScalarField K2("K2", 0.5*magSqr(U2)); + volScalarField K1("K" + phase1Name, 0.5*magSqr(U1)); + volScalarField K2("K" + phase2Name, 0.5*magSqr(U2)); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H index 8f05a1d1b3..0a782ef51e 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H @@ -151,7 +151,7 @@ ( IOobject ( - "nut2", + "nut" + phase2Name, runTime.timeName(), mesh, IOobject::NO_READ, @@ -165,13 +165,13 @@ ( IOobject ( - "nuEff1", + "nuEff" + phase1Name, runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), - sqr(Ct)*nut2 + nu1 + sqr(Ct)*nut2 + thermo1.mu()/rho1 ); Info<< "Calculating field nuEff2\n" << endl; @@ -179,11 +179,11 @@ ( IOobject ( - "nuEff2", + "nuEff" + phase2Name, runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), - nut2 + nu2 + nut2 + thermo2.mu()/rho2 ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H index d6ccf90289..d53bec5ea4 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H @@ -38,12 +38,12 @@ volScalarField heatTransferCoeff volVectorField Ur(U1 - U2); volScalarField magUr(mag(Ur) + residualSlip); - if (dispersedPhase == "1") + if (dispersedPhase == phase1Name) { dragCoeff = drag1->K(magUr); heatTransferCoeff = heatTransfer1->K(magUr); } - else if (dispersedPhase == "2") + else if (dispersedPhase == phase2Name) { dragCoeff = drag2->K(magUr); heatTransferCoeff = heatTransfer2->K(magUr); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options index 92fe0b0d74..f031e05898 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I../phaseModel/lnInclude LIB_LIBS = \ - -lcompressiblePhaseModel - + -lcompressiblePhaseModel \ + -lfluidThermophysicalModels \ + -lspecie diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C index 9bc6cadaf6..95ade6ba05 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,8 +72,7 @@ Foam::tmp Foam::heatTransferModels::RanzMarshall::K ) const { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - dimensionedScalar Prb = - phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa(); + volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa()); volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb)); return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d()); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options index 2fcce9913d..7fdabf5730 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options @@ -1,5 +1,6 @@ EXE_INC = \ -I$(LIB_SRC)/foam/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I../phaseModel/lnInclude \ -I../interfacialModels/lnInclude diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index 7d0767914e..18db3a70df 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel phi1_(phase1.phi()), draga_(draga), - rho1_(phase1.rho()), - nu1_(phase1.nu()), + rho1_(phase1.rho()[0]), //***HGW + nu1_(phase1.nu()()[0]), //***HGW kineticTheoryProperties_ ( @@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel ( IOobject ( - "mu1", + "mu" + phase1.name(), U1_.time().timeName(), U1_.mesh(), IOobject::NO_READ, diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index f7010b1c45..f9a7807d00 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -1,6 +1,6 @@ { - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; + rho1 = thermo1.rho(); + rho2 = thermo2.rho(); surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField alpha2f(scalar(1) - alpha1f); @@ -11,10 +11,10 @@ surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1)); surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2)); - volVectorField HbyA1("HbyA1", U1); + volVectorField HbyA1("HbyA" + phase1Name, U1); HbyA1 = rAU1*U1Eqn.H(); - volVectorField HbyA2("HbyA2", U2); + volVectorField HbyA2("HbyA" + phase2Name, U2); HbyA2 = rAU2*U2Eqn.H(); mrfZones.absoluteFlux(phi1.oldTime()); @@ -38,14 +38,14 @@ surfaceScalarField phiHbyA1 ( - "phiHbyA1", + "phiHbyA" + phase1Name, (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1) ); surfaceScalarField phiHbyA2 ( - "phiHbyA2", + "phiHbyA" + phase2Name, (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2) ); @@ -96,8 +96,16 @@ //} //else { - surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1); - surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2); + surfaceScalarField phid1 + ( + "phid" + phase1Name, + fvc::interpolate(psi1)*phi1 + ); + surfaceScalarField phid2 + ( + "phid" + phase2Name, + fvc::interpolate(psi2)*phi2 + ); pEqnComp1 = fvc::ddt(rho1) + psi1*correction(fvm::ddt(p)) @@ -174,13 +182,15 @@ p = max(p, pMin); - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; + thermo1.correct(); + thermo2.correct(); + rho1 = thermo1.rho(); + rho2 = thermo2.rho(); K1 = 0.5*magSqr(U1); K2 = 0.5*magSqr(U2); - //***HGW if (thermo.dpdt()) + if (thermo1.dpdt()) { dpdt = fvc::ddt(p); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options index 0ec1139209..e441b0417b 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options @@ -1,6 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude LIB_LIBS = \ - -lincompressibleTransportModels + -lincompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C index d80758ac8c..64d59b436d 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel const word& phaseName ) : - dict_ + volScalarField ( - transportProperties.subDict("phase" + phaseName) + IOobject + ( + "alpha" + phaseName, + mesh.time().timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("alpha", dimless, 0) ), name_(phaseName), - nu_ - ( - "nu", - dimensionSet(0, 2, -1, 0, 0), - dict_.lookup("nu") - ), - kappa_ - ( - "kappa", - dimensionSet(1, 1, -3, -1, 0), - dict_.lookup("kappa") - ), - Cp_ - ( - "Cp", - dimensionSet(0, 2, -2, -1, 0), - dict_.lookup("Cp") - ), - rho_ - ( - "rho", - dimDensity, - dict_.lookup("rho") - ), + phaseDict_(transportProperties.subDict(phaseName)), + thermo_(rhoThermo::New(mesh, phaseName)), U_ ( IOobject @@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel mesh ) { + thermo_->validate("phaseModel " + phaseName, "h", "e"); + const word phiName = "phi" + phaseName; IOobject phiHeader @@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel dPtr_ = diameterModel::New ( - dict_, + phaseDict_, *this ); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H index 90b2d2184e..18f5a47e46 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,6 +36,7 @@ SourceFiles #include "dimensionedScalar.H" #include "volFields.H" #include "surfaceFields.H" +#include "rhoThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -50,25 +51,18 @@ class diameterModel; \*---------------------------------------------------------------------------*/ class phaseModel +: + public volScalarField { // Private data - dictionary dict_; - //- Name of phase word name_; - //- Kinematic viscosity - dimensionedScalar nu_; + dictionary phaseDict_; - //- Thermal conductivity - dimensionedScalar kappa_; - - //- Heat capacity - dimensionedScalar Cp_; - - //- Density - dimensionedScalar rho_; + //- Thermophysical properties + autoPtr thermo_; //- Velocity volVectorField U_; @@ -116,24 +110,34 @@ public: tmp d() const; - const dimensionedScalar& nu() const + tmp nu() const { - return nu_; + return thermo_->mu()/thermo_->rho(); } - const dimensionedScalar& kappa() const + tmp kappa() const { - return kappa_; + return thermo_->kappa(); } - const dimensionedScalar& Cp() const + tmp Cp() const { - return Cp_; + return thermo_->Cp(); } - const dimensionedScalar& rho() const + const volScalarField& rho() const { - return rho_; + return thermo_->rho(); + } + + const rhoThermo& thermo() const + { + return thermo_(); + } + + rhoThermo& thermo() + { + return thermo_(); } const volVectorField& U() const diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H index 7f53411529..bc9a07b0a8 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H @@ -61,4 +61,4 @@ if (turbulence) #include "wallViscosity.H" } -nuEff2 = nut2 + nu2; +nuEff2 = nut2 + thermo2.mu()/rho2; diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H index c91cce9a00..d85181cba2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H @@ -4,7 +4,6 @@ scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar Cmu75 = ::pow(Cmu.value(), 0.75); scalar kappa_ = kappa.value(); - scalar nu2_ = nu2.value(); const fvPatchList& patches = mesh.boundary(); @@ -30,6 +29,8 @@ forAll(patches, patchi) { const fvPatch& currPatch = patches[patchi]; + const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi]; + const scalarField& rho2_ = rho2.boundaryField()[patchi]; if (isA(currPatch)) { @@ -52,7 +53,8 @@ /(kappa_*y[patchi][facei]); G[faceCelli] += - (nut2w[facei] + nu2_)*magFaceGradU[facei] + (nut2w[facei] + mu2_[facei]/rho2_[facei]) + *magFaceGradU[facei] *Cmu25*::sqrt(k[faceCelli]) /(kappa_*y[patchi][facei]); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H index b153a36014..9aa032149c 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H @@ -2,13 +2,14 @@ scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar kappa_ = kappa.value(); scalar E_ = E.value(); - scalar nu2_ = nu2.value(); const fvPatchList& patches = mesh.boundary(); forAll(patches, patchi) { const fvPatch& currPatch = patches[patchi]; + const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi]; + const scalarField& rho2_ = rho2.boundaryField()[patchi]; if (isA(currPatch)) { @@ -20,11 +21,14 @@ // calculate yPlus scalar yPlus = - Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_; + Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) + /(mu2_[facei]/rho2_[facei]); if (yPlus > 11.6) { - nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1); + nutw[facei] = + (mu2_[facei]/rho2_[facei]) + *(yPlus*kappa_/::log(E_*yPlus) - 1); } else { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H index 34fdc24afc..ddb923cf87 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H @@ -13,7 +13,7 @@ surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf()); phir += phipp; - phic += fvc::interpolate(alpha1)*phipp; + phic += alpha1f*phipp; } for (int acorr=0; acorr 0) { - ppMagf = rAU1f*fvc::interpolate - ( - (1.0/(rho1*(alpha1 + scalar(0.0001)))) - *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) - ); + surfaceScalarField alpha1f(fvc::interpolate(alpha1)); + + // ppMagf = rAU1f*fvc::interpolate + // ( + // (1.0/(rho1*(alpha1 + scalar(0.0001)))) + // *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) + // ); + ppMagf = + rAU1f/(alpha1f + scalar(0.0001)) + *(g0/rho1)*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax); fvScalarMatrix alpha1Eqn ( fvm::ddt(alpha1) - fvc::ddt(alpha1) - fvm::laplacian ( - (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, + alpha1f*ppMagf, alpha1, "laplacian(alpha1PpMag,alpha1)" ) diff --git a/applications/test/Circulator/Test-Circulator.C b/applications/test/Circulator/Test-Circulator.C index 19a21f6b79..2d8ecd3765 100644 --- a/applications/test/Circulator/Test-Circulator.C +++ b/applications/test/Circulator/Test-Circulator.C @@ -38,94 +38,6 @@ Description using namespace Foam; -// return -// 0: no match -// +1: identical -// -1: same face, but different orientation -label compare(const face& a, const face& b) -{ - // Basic rule: we assume that the sequence of labels in each list - // will be circular in the same order (but not necessarily in the - // same direction or from the same starting point). - - // Trivial reject: faces are different size - label sizeA = a.size(); - label sizeB = b.size(); - - if (sizeA != sizeB || sizeA == 0) - { - return 0; - } - - const_circulator aCirc(a); - const_circulator bCirc(b); - - // Rotate face b until its element matches the starting element of face a. - do - { - if (aCirc() == bCirc()) - { - // Set bCirc fulcrum to its iterator and increment the iterators - bCirc.setFulcrumToIterator(); - ++aCirc; - ++bCirc; - - break; - } - } while (bCirc.circulate(CirculatorBase::CLOCKWISE)); - - // Look forwards around the faces for a match - do - { - if (aCirc() != bCirc()) - { - break; - } - } - while - ( - aCirc.circulate(CirculatorBase::CLOCKWISE), - bCirc.circulate(CirculatorBase::CLOCKWISE) - ); - - // If the circulator has stopped then faces a and b matched. - if (!aCirc.circulate()) - { - return 1; - } - else - { - // Reset the circulators back to their fulcrum - aCirc.setIteratorToFulcrum(); - bCirc.setIteratorToFulcrum(); - ++aCirc; - --bCirc; - } - - // Look backwards around the faces for a match - do - { - if (aCirc() != bCirc()) - { - break; - } - } - while - ( - aCirc.circulate(CirculatorBase::CLOCKWISE), - bCirc.circulate(CirculatorBase::ANTICLOCKWISE) - ); - - // If the circulator has stopped then faces a and b matched. - if (!aCirc.circulate()) - { - return -1; - } - - return 0; -} - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: @@ -184,40 +96,44 @@ int main(int argc, char *argv[]) Info<< nl << nl << "Compare two faces: " << endl; face a(identity(5)); - Info<< "Compare " << a << " and " << a << " Match = " << compare(a, a) + Info<< "Compare " << a << " and " << a << " Match = " << face::compare(a, a) << endl; face b(reverseList(a)); - Info<< "Compare " << a << " and " << b << " Match = " << compare(a, b) + Info<< "Compare " << a << " and " << b << " Match = " << face::compare(a, b) << endl; face c(a); c[4] = 3; - Info<< "Compare " << a << " and " << c << " Match = " << compare(a, c) + Info<< "Compare " << a << " and " << c << " Match = " << face::compare(a, c) << endl; face d(rotateList(a, 2)); - Info<< "Compare " << a << " and " << d << " Match = " << compare(a, d) + Info<< "Compare " << a << " and " << d << " Match = " << face::compare(a, d) << endl; face g(labelList(5, 1)); face h(g); - Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h) + Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h) << endl; g[0] = 2; h[3] = 2; - Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h) + Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h) << endl; g[4] = 3; h[4] = 3; - Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h) + Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h) << endl; face face1(identity(1)); Info<< "Compare " << face1 << " and " << face1 - << " Match = " << compare(face1, face1) << endl; + << " Match = " << face::compare(face1, face1) << endl; + + face face2(identity(1)+1); + Info<< "Compare " << face1 << " and " << face2 + << " Match = " << face::compare(face1, face2) << endl; Info<< nl << nl << "Zero face" << nl << endl; diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict index 60f4cb2f43..340c33049f 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseDict +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict @@ -17,7 +17,7 @@ FoamFile collapseEdgesCoeffs { // Edges shorter than this absolute value will be merged - minimumEdgeLength 1e-8; + minimumEdgeLength 1e-6; // The maximum angle between two edges that share a point attached to // no other edges @@ -25,7 +25,7 @@ collapseEdgesCoeffs // The amount that minimumEdgeLength will be reduced by for each // edge if that edge's collapse generates a poor quality face - reductionFactor 0.5; + reductionFactor 0.5; } @@ -67,14 +67,17 @@ meshQualityCoeffs { // Name of the dictionary that has the mesh quality coefficients used // by motionSmoother::checkMesh - meshQualityCoeffDict meshQualityDict; + #include "meshQualityDict"; + + // Maximum number of smoothing iterations for the reductionFactors + maximumSmoothingIterations 2; // Maximum number of outer iterations is mesh quality checking is enabled - maximumIterations 10; + maximumIterations 10; // Maximum number of iterations deletion of a point can cause a bad face // to be constructed before it is forced to not be deleted - maxPointErrorCount 5; + maxPointErrorCount 5; } diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict index ffa2db8b19..6587c84298 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict +++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMeshDict @@ -70,6 +70,9 @@ extrudeModel linearNormal; //- Extrudes into sphere around (0 0 0) // extrudeModel linearRadial; +//- Extrudes into sphere around (0 0 0) with specified radii +//extrudeModel radial; + //- Extrudes into sphere with grading according to pressure (atmospherics) // extrudeModel sigmaRadial; @@ -98,6 +101,14 @@ linearDirectionCoeffs linearRadialCoeffs { R 0.1; + // Optional inner radius + Rsurface 0.01; +} + +radialCoeffs +{ + // Radii specified through interpolation table + R table ((0 0.01)(3 0.03)(10 0.1)); } sigmaRadialCoeffs @@ -107,4 +118,5 @@ sigmaRadialCoeffs pStrat 1; } + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options index f1db60d53f..fe1283a73a 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options @@ -1,5 +1,4 @@ EXE_INC = \ - /* -DFULLDEBUG -g -O0 */ \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C index 7979fa966a..1b83599b31 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C @@ -26,12 +26,6 @@ License #include "patchToPoly2DMesh.H" #include "PatchTools.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // - - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::patchToPoly2DMesh::flipFaceOrder() @@ -275,12 +269,9 @@ void Foam::patchToPoly2DMesh::createPolyMeshComponents() addPatchFacesToOwner(); } -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -//- Construct from a primitivePatch Foam::patchToPoly2DMesh::patchToPoly2DMesh ( const MeshedSurface& patch, @@ -321,13 +312,5 @@ void Foam::patchToPoly2DMesh::createMesh() } } -// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * // - // ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H index ea27dc6f53..b7c81e68a2 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H @@ -51,6 +51,7 @@ class patchToPoly2DMesh { // Private data + // Reference to the meshed surface const MeshedSurface& patch_; const wordList& patchNames_; @@ -62,11 +63,16 @@ class patchToPoly2DMesh const EdgeMap