From 63b2aa37ead6cfc6c9553b633cf78bee956dc82d Mon Sep 17 00:00:00 2001 From: danielque Date: Thu, 12 Dec 2019 12:41:51 +0100 Subject: [PATCH 01/75] start implementation of recurrence chemistry solver not quite sure yet what to do with rhoEqn --- .../rcfdemSolverRhoSteadyPimpleChem/EEqn.H | 55 +++ .../Make/files | 3 + .../Make/options | 52 +++ .../rcfdemSolverRhoSteadyPimpleChem/UEqn.H | 33 ++ .../rcfdemSolverRhoSteadyPimpleChem/YEqn.H | 81 ++++ .../createFieldRefs.H | 2 + .../createFields.H | 416 ++++++++++++++++++ .../rcfdemSolverRhoSteadyPimpleChem/limitP.H | 2 + .../rcfdemSolverRhoSteadyPimpleChem/limitU.H | 11 + .../rcfdemSolverRhoSteadyPimpleChem/molConc.H | 12 + .../monitorMass.H | 7 + .../rcfdemSolverRhoSteadyPimpleChem/pEqn.H | 91 ++++ .../rcfdemSolverRhoSteadyPimpleChem.C | 218 +++++++++ .../rcfdemSolverRhoSteadyPimpleChem/rhoEqn.H | 21 + .../updateFields.H | 8 + etc/solver-list.txt | 2 + 16 files changed, 1014 insertions(+) create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rhoEqn.H create mode 100644 applications/solvers/rcfdemSolverRhoSteadyPimpleChem/updateFields.H diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H new file mode 100644 index 00000000..692f4873 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H @@ -0,0 +1,55 @@ +// contributions to internal energy equation can be found in +// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998 +{ + // dim he = J / kg + volScalarField& he = thermo.he(); + particleCloud.energyContributions(Qsource); + particleCloud.energyCoefficients(QCoeff); + + addSource = + ( + he.name() == "e" + ? + fvc::div(phi, K) + + fvc::div + ( + fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U), + p, + "div(phiv,p)" + ) + : fvc::div(phi, K) + ); + + Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); + + + fvScalarMatrix EEqn + ( + fvm::div(phi, he) + + addSource + // net heat transfer from particles to fluid + - Qsource + - fvm::Sp(QCoeff/Cpv, he) + // - fvm::laplacian(voidfractionRec*kf/Cpv,he) + - fvm::laplacian(voidfractionRec*thCond/Cpv,he) + == + fvOptions(rho, he) + ); + + EEqn.relax(); + + fvOptions.constrain(EEqn); + + EEqn.solve(); + + fvOptions.correct(he); + + thermo.correct(); + + Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl; + + + particleCloud.clockM().start(31,"energySolve"); + particleCloud.solve(); + particleCloud.clockM().stop("energySolve"); +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files new file mode 100644 index 00000000..a1096b62 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files @@ -0,0 +1,3 @@ +rcfdemSolverRhoSteadyPimpleChem.C + +EXE=$(CFDEM_APP_DIR)/rcfdemSolverRhoSteadyPimpleChem diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options new file mode 100644 index 00000000..9290ba36 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options @@ -0,0 +1,52 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION))) +PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR) +PFLAGS+= -Dcompre + +EXE_INC = \ + $(PFLAGS) \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/combustionModels/lnInclude \ + -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lrecurrence \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lfvOptions \ + -l$(CFDEM_LIB_COMP_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) \ + -lreactionThermophysicalModels \ + -lchemistryModel \ + -lradiationModels \ + -lregionModels \ + -lsurfaceFilmModels \ + -lODE \ + -lcombustionModels diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H new file mode 100644 index 00000000..9e75193e --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H @@ -0,0 +1,33 @@ +// Solve the Momentum equation +particleCloud.otherForces(fOther); + +tmp tUEqn +( + fvm::div(phi, U) + + particleCloud.divVoidfractionTau(U, voidfractionRec) + + fvm::Sp(Ksl,U) + - fOther + == + fvOptions(rho, U) +); +fvVectorMatrix& UEqn = tUEqn.ref(); + +UEqn.relax(); + +fvOptions.constrain(UEqn); + +if (modelType=="B" || modelType=="Bfull") +{ + solve(UEqn == -fvc::grad(p)+ Ksl*UsRec); +} +else +{ + solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec); +} + +//U.relax(); +#include "limitU.H" + +fvOptions.correct(U); + +K = 0.5*magSqr(U); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H new file mode 100644 index 00000000..aaee90a2 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H @@ -0,0 +1,81 @@ +particleCloud.clockM().start(29,"Y"); + +tmp > mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); + +{ + combustion->correct(); +#if OPENFOAM_VERSION_MAJOR < 5 + dQ = combustion->dQ(); +#else + Qdot = combustion->Qdot(); +#endif + label inertIndex = -1; + volScalarField Yt(0.0*Y[0]); + + forAll(Y, i) + { + if (Y[i].name() == inertSpecie) inertIndex = i; + if (Y[i].name() != inertSpecie || propagateInertSpecie) + { + volScalarField& Yi = Y[i]; + + fvScalarMatrix YiEqn + ( + //fvm::ddt(rhoeps, Yi) + //+ + mvConvection->fvmDiv(phi, Yi) + - fvm::laplacian(voidfraction*turbulence->muEff(), Yi) + == + combustion->R(Yi) + + particleCloud.chemistryM(0).Smi(i) + + fvOptions(rho, Yi) + ); + + YiEqn.relax(); + + fvOptions.constrain(YiEqn); + + YiEqn.solve(mesh.solver("Yi")); + + fvOptions.correct(Yi); + + Yi.max(0.0); + if (Y[i].name() != inertSpecie) Yt += Yi; + } + } + + if (inertIndex!=-1) + { + Y[inertIndex].max(inertLowerBound); + Y[inertIndex].min(inertUpperBound); + } + + if (propagateInertSpecie) + { + if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL); + forAll(Y,i) + { + if (i!=inertIndex) + { + volScalarField& Yi = Y[i]; + Yi = Yi/(Yt+VSMALL); + } + } + } + else + { + Y[inertIndex] = scalar(1) - Yt; + Y[inertIndex].max(0.0); + } +} + +particleCloud.clockM().stop("Y"); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H new file mode 100644 index 00000000..5842906a --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H @@ -0,0 +1,2 @@ +const volScalarField& T = thermo.T(); +const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H new file mode 100644 index 00000000..23cb60fc --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H @@ -0,0 +1,416 @@ +Info<< "Reading thermophysical properties\n" << endl; + +#if OPENFOAM_VERSION_MAJOR < 6 + Info<< "Creating combustion model\n" << endl; + autoPtr combustion + ( + combustionModels::rhoCombustionModel::New(mesh) + ); + rhoReactionThermo& thermo = combustion->thermo(); +#else + Info<< "Reading thermophysical properties\n" << endl; + autoPtr pThermo(rhoReactionThermo::New(mesh)); + rhoReactionThermo& thermo = pThermo(); +#endif + thermo.validate(args.executable(), "h", "e"); + + basicSpecieMixture& composition = thermo.composition(); + PtrList& Y = composition.Y(); + + // read molecular weight +#if OPENFOAM_VERSION_MAJOR < 6 + volScalarField W(composition.W()); +#else + volScalarField W(thermo.W()); +#endif + + bool propagateInertSpecie = true; + + const word inertSpecie(thermo.lookup("inertSpecie")); + + const scalar inertLowerBound(thermo.lookupOrDefault("inertLowerBound",0.0)); + + const scalar inertUpperBound(thermo.lookupOrDefault("inertUpperBound",1.0)); + + 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(); + + multivariateSurfaceInterpolationScheme::fieldTable fields; + + forAll(Y, i) + { + fields.add(Y[i]); + } + fields.add(thermo.he()); + + Info<< "Reading field rho\n" << endl; + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl; + volScalarField voidfraction + ( + IOobject + ( + "voidfraction", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + volScalarField voidfractionRec + ( + IOobject + ( + "voidfractionRec", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + voidfraction + ); + + volScalarField addSource + ( + IOobject + ( + "addSource", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0) + ); + + Info<< "\nCreating fluid-particle heat flux field\n" << endl; + volScalarField Qsource + ( + IOobject + ( + "Qsource", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0) + ); + + Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl; + volScalarField QCoeff + ( + IOobject + ( + "QCoeff", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0) + ); + + Info<< "\nCreating thermal conductivity field\n" << endl; + volScalarField thCond + ( + IOobject + ( + "thCond", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0), + "zeroGradient" + ); + + Info<< "\nCreating heat capacity field\n" << endl; + volScalarField Cpv + ( + IOobject + ( + "Cpv", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0) + ); + + Info<< "\nCreating body force field\n" << endl; + volVectorField fOther + ( + IOobject + ( + "fOther", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) + ); + + Info<< "Reading/calculating face flux field phi\n" << endl; + surfaceScalarField phi + ( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(rho*U*voidfraction) & mesh.Sf() + ); + + dimensionedScalar rhoMax + ( + dimensionedScalar::lookupOrDefault + ( + "rhoMax", + pimple.dict(), + dimDensity, + GREAT + ) + ); + + dimensionedScalar rhoMin + ( + dimensionedScalar::lookupOrDefault + ( + "rhoMin", + pimple.dict(), + dimDensity, + 0 + ) + ); + + dimensionedScalar pMax + ( + dimensionedScalar::lookupOrDefault + ( + "pMax", + pimple.dict(), + dimPressure, + GREAT + ) + ); + + dimensionedScalar pMin + ( + dimensionedScalar::lookupOrDefault + ( + "pMin", + pimple.dict(), + dimPressure, + -GREAT + ) + ); + + dimensionedScalar UMax + ( + dimensionedScalar::lookupOrDefault + ( + "UMax", + pimple.dict(), + dimVelocity, + -1.0 + ) + ); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); + +#if OPENFOAM_VERSION_MAJOR >= 6 + Info<< "Creating combustion model\n" << endl; + autoPtr> combustion + ( + CombustionModel::New(thermo, turbulence()) + ); +#endif + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, pimple.dict(), pRefCell, pRefValue); + + mesh.setFluxRequired(p.name()); + + Info<< "Creating field dpdt\n" << endl; + volScalarField dpdt + ( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("dpdt", p.dimensions()/dimTime, 0) + ); + + Info<< "Creating field kinetic energy K\n" << endl; + volScalarField K("K", 0.5*magSqr(U)); + +#if OPENFOAM_VERSION_MAJOR < 5 + volScalarField dQ + ( + IOobject + ( + "dQ", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) + ); +#else + volScalarField Qdot + ( + IOobject + ( + "Qdot", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0) + ); +#endif + + Info<< "\nReading momentum exchange field Ksl\n" << endl; + volScalarField Ksl + ( + IOobject + ( + "Ksl", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0) + ); + + + Info<< "Reading particle velocity field Us\n" << endl; + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + volScalarField molarConc + ( + IOobject + ( + "molarConc", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0) + ); + + volScalarField dSauter + ( + IOobject + ( + "dSauter", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0) + ); + + volVectorField UsRec + ( + IOobject + ( + "UsRec", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + Us + ); + + + dimensionedScalar kf("0", dimensionSet(1, 1, -3, -1, 0, 0, 0), 0.026); + +//=============================== diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H new file mode 100644 index 00000000..f0971a14 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H @@ -0,0 +1,2 @@ +p = max(p, pMin); +p = min(p, pMax); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H new file mode 100644 index 00000000..f699b25c --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H @@ -0,0 +1,11 @@ +if (UMax.value() > 0) +{ + forAll(U,cellI) + { + scalar mU(mag(U[cellI])); + if (mU > UMax.value()) + { + U[cellI] *= UMax.value() / mU; + } + } +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H new file mode 100644 index 00000000..dc70981e --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H @@ -0,0 +1,12 @@ +{ + molarConc = 0.0 * molarConc; + forAll(Y, i) + { + volScalarField& Yi = Y[i]; + dimensionedScalar mi("mi",dimensionSet(1, 0, 0, 0, -1),composition.W(i)); + mi /= 1000.0; // g to kg + molarConc += rho * Yi / mi; + } +} + +// ************************************************************************* // diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H new file mode 100644 index 00000000..5940864f --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H @@ -0,0 +1,7 @@ +{ + m=gSum(rhoeps*1.0*rhoeps.mesh().V()); + if(counter==0) m0=m; + counter++; + Info << "\ncurrent gas mass = " << m << "\n" << endl; + Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl; +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H new file mode 100644 index 00000000..6b40829c --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H @@ -0,0 +1,91 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU)); +if (modelType=="A") +{ + rhorAUf *= fvc::interpolate(voidfractionRec); +} +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf()); + + +if (pimple.transonic()) +{ +// transonic version not implemented yet +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + fvc::flux(rhoeps*HbyA) + ) + ); + + // flux without pressure gradient contribution + phi = phiHbyA + phiUs; + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rhoeps, U, phi, rhorAUf); + + volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p); + + while (pimple.correctNonOrthogonal()) + { + // Pressure corrector + fvScalarMatrix pEqn + ( + fvc::div(phi) + - fvm::laplacian(rhorAUf, p) + == + fvm::Sp(SmbyP, p) + + fvOptions(psi, p, rho.name()) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi += pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrsPU.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +#include "limitP.H" + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); +Info<< "rho max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + +if (modelType=="A") +{ + U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec); +} +else +{ + U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec); +} + +#include "limitU.H" + +U.correctBoundaryConditions(); +fvOptions.correct(U); +K = 0.5*magSqr(U); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C new file mode 100644 index 00000000..9e22b7f8 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C @@ -0,0 +1,218 @@ +/*---------------------------------------------------------------------------*\ +License + + This 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. + + This code 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 this code. If not, see . + + Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria + +Application + rcfdemSolverRhoSteadyPimpleChem + +Description + Transient (DEM) + steady-state (CFD) solver for compressible flow using the + flexible PIMPLE (PISO-SIMPLE) algorithm. Particle-motion is obtained from + a recurrence process. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x, + where additional functionality for CFD-DEM coupling is added. +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +//#include "psiThermo.H" +#include "turbulentFluidThermoModel.H" +#if OPENFOAM_VERSION_MAJOR < 6 +#include "rhoCombustionModel.H" +#else +#include "rhoReactionThermo.H" +#include "CombustionModel.H" +#endif +#include "bound.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" + +#include "cfdemCloudRec.H" +#include "recBase.H" +#include "recModel.H" +#include "recPath.H" + +#include "cfdemCloudEnergy.H" +#include "implicitCouple.H" +#include "clockModel.H" +#include "smoothingModel.H" +#include "forceModel.H" +#include "thermCondModel.H" +#include "energyModel.H" +#include "chemistryModel.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 "createRDeltaT.H" + + #include "initContinuityErrs.H" + #include "createFields.H" + //#include "createFieldRefs.H" + #include "createFvOptions.H" + + // create cfdemCloud + //#include "readGravitationalAcceleration.H" + cfdemCloudRec particleCloud(mesh); + #include "checkModelType.H" + recBase recurrenceBase(mesh); + #include "updateFields.H" + + turbulence->validate(); + //#include "compressibleCourantNo.H" + //#include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + label recTimeIndex = 0; + scalar recTimeStep = recurrenceBase.recM().recTimeStep(); + scalar startTime = runTime.startTime().value(); + + const IOdictionary& couplingProps = particleCloud.couplingProperties(); + label nEveryFlow(couplingProps.lookupOrDefault