diff --git a/applications/solvers/cfdemSolverPiso/UEqn.H b/applications/solvers/cfdemSolverPiso/UEqn.H index 2f26ac91..21165a36 100644 --- a/applications/solvers/cfdemSolverPiso/UEqn.H +++ b/applications/solvers/cfdemSolverPiso/UEqn.H @@ -1,8 +1,11 @@ +particleCloud.otherForces(fOther); + fvVectorMatrix UEqn ( fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U) + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U) + particleCloud.divVoidfractionTau(U, voidfraction) + - fOther/rho == fvOptions(U) - fvm::Sp(Ksl/rho,U) diff --git a/applications/solvers/cfdemSolverPiso/createFields.H b/applications/solvers/cfdemSolverPiso/createFields.H index 36246301..72f12354 100644 --- a/applications/solvers/cfdemSolverPiso/createFields.H +++ b/applications/solvers/cfdemSolverPiso/createFields.H @@ -46,6 +46,21 @@ //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.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<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl; volScalarField voidfraction ( diff --git a/applications/solvers/cfdemSolverPisoFreeStreaming/Make/files b/applications/solvers/cfdemSolverPisoFreeStreaming/Make/files new file mode 100644 index 00000000..a91d0a6c --- /dev/null +++ b/applications/solvers/cfdemSolverPisoFreeStreaming/Make/files @@ -0,0 +1,3 @@ +cfdemSolverPisoFreeStreaming.C + +EXE=$(CFDEM_APP_DIR)/cfdemSolverPisoFreeStreaming diff --git a/applications/solvers/cfdemSolverPisoFreeStreaming/Make/options b/applications/solvers/cfdemSolverPisoFreeStreaming/Make/options new file mode 100644 index 00000000..ef370ae5 --- /dev/null +++ b/applications/solvers/cfdemSolverPisoFreeStreaming/Make/options @@ -0,0 +1,28 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +EXE_INC = \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(FOAM_SOLVERS)/incompressible/pisoFoam \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools \ + -lfvOptions \ + -lsampling \ + -l$(CFDEM_LIB_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/cfdemSolverPisoFreeStreaming/cfdemSolverPisoFreeStreaming.C b/applications/solvers/cfdemSolverPisoFreeStreaming/cfdemSolverPisoFreeStreaming.C new file mode 100644 index 00000000..7e370549 --- /dev/null +++ b/applications/solvers/cfdemSolverPisoFreeStreaming/cfdemSolverPisoFreeStreaming.C @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright (C) 1991-2009 OpenCFD Ltd. + Copyright (C) 2009-2012 JKU, Linz + Copyright (C) 2012- DCS Computing GmbH,Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling 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. + + CFDEMcoupling 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 CFDEMcoupling. If not, see . + +Application + cfdemSolverPisoFreeStreaming + +Description + Transient solver for incompressible flow. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, + where additional functionality for CFD-DEM coupling is added. + the particles follow the fluid velocity +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "pisoControl.H" +#include "fvOptions.H" + +#include "cfdemCloudRec.H" + +#include "cfdemCloud.H" +#include "implicitCouple.H" +#include "clockModel.H" +#include "smoothingModel.H" +#include "forceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + + cfdemCloudRec particleCloud(mesh); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + Info<< "\nStarting time loop\n" << endl; + while (runTime.loop()) + { + + particleCloud.clockM().start(1,"Global"); + + Info<< "Time = " << runTime.timeName() << nl << endl; + + #include "CourantNo.H" + + // do particle stuff + particleCloud.clockM().start(2,"Coupling"); + + particleCloud.evolve(voidfraction,Us,U); + + + particleCloud.clockM().stop("Coupling"); + + particleCloud.clockM().start(26,"Flow"); + + if(particleCloud.solveFlow()) + { + // Pressure-velocity PISO corrector + { + // Momentum predictor + #include "UEqn.H" + + // --- PISO loop + + while (piso.correct()) + { + #include "pEqn.H" + } + } + + laminarTransport.correct(); + turbulence->correct(); + } + else + { + Info << "skipping flow solution." << endl; + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + particleCloud.clockM().stop("Flow"); + particleCloud.clockM().stop("Global"); + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverPisoFreeStreaming/createFields.H b/applications/solvers/cfdemSolverPisoFreeStreaming/createFields.H new file mode 100644 index 00000000..598446e0 --- /dev/null +++ b/applications/solvers/cfdemSolverPisoFreeStreaming/createFields.H @@ -0,0 +1,110 @@ + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading physical velocity field U" << endl; + Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//=============================== +// particle interaction modelling +//=============================== + + Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl; + volScalarField voidfraction + ( + IOobject + ( + "voidfraction", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "\nCreating density field rho\n" << endl; + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading particle velocity field Us\n" << endl; + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//=============================== + +//# include "createPhi.H" +#ifndef createPhi_H +#define createPhi_H +Info<< "Reading/calculating face flux field phi\n" << endl; +surfaceScalarField phi +( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(U) & mesh.Sf() +); +#endif + + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); + + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + +#include "createMRF.H" diff --git a/applications/solvers/cfdemSolverPisoScalar/createFields.H b/applications/solvers/cfdemSolverPisoScalar/createFields.H index f7779bec..9cd6599f 100644 --- a/applications/solvers/cfdemSolverPisoScalar/createFields.H +++ b/applications/solvers/cfdemSolverPisoScalar/createFields.H @@ -46,6 +46,21 @@ //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.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<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl; volScalarField voidfraction ( diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index 022e933d..afce5d6d 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H @@ -20,21 +20,23 @@ Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); - // correct source for the thermodynamic reference temperature - dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL)); - Qsource += QCoeff*Tref; +// For implict T terms in the energy/enthalpy transport equation, use +// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1. +// This formula is valid for ideal gases with e=e(T) and h=h(T). For +// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction +// terms accounting for pressure variations. fvScalarMatrix EEqn ( fvm::ddt(rhoeps, he) + fvm::div(phi, he) + addSource - // net heat transfer from particles to fluid - Qsource + - QCoeff*T - fvm::Sp(QCoeff/Cpv, he) - // thermal conduction of the fluid with effective conductivity + + QCoeff/Cpv*he + - fvc::laplacian(voidfraction*thCond,T) - fvm::laplacian(voidfraction*thCond/Cpv,he) - // + particle-fluid energy transfer due to work - // + fluid energy dissipation due to shearing + + fvc::laplacian(voidfraction*thCond/Cpv,he) == fvOptions(rho, he) ); diff --git a/applications/solvers/rcfdemSolverBase/createFields.H b/applications/solvers/rcfdemSolverBase/createFields.H index 2d002e5e..7088bab5 100644 --- a/applications/solvers/rcfdemSolverBase/createFields.H +++ b/applications/solvers/rcfdemSolverBase/createFields.H @@ -37,12 +37,20 @@ "URec", runTime.timeName(), mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE ), - mesh + mesh, + dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero) ); + bool updateURec = false; + if (URec.headerOk()) + { + updateURec = true; + URec.writeOpt() = IOobject::AUTO_WRITE; + } + volScalarField voidfractionRec ( IOobject @@ -50,12 +58,20 @@ "voidfractionRec", runTime.timeName(), mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE ), - mesh + mesh, + dimensionedScalar("voidfractionRec", dimensionSet(0, 0, 0, 0, 0), 1.0) ); + bool updateVoidfractionRec = false; + if (voidfractionRec.headerOk()) + { + updateVoidfractionRec = true; + voidfractionRec.writeOpt() = IOobject::AUTO_WRITE; + } + volVectorField UsRec ( IOobject @@ -63,12 +79,20 @@ "UsRec", runTime.timeName(), mesh, - IOobject::MUST_READ, + IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - mesh + mesh, + dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero) ); + bool updateUsRec = false; + if (UsRec.headerOk()) + { + updateUsRec = true; + UsRec.writeOpt() = IOobject::AUTO_WRITE; + } + // calculated fields Info << "\nCreating fields subject to calculation\n" << endl; volScalarField voidfraction @@ -78,7 +102,7 @@ "voidfraction", runTime.timeName(), mesh, - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), voidfractionRec @@ -91,7 +115,7 @@ "Us", runTime.timeName(), mesh, - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), UsRec @@ -111,11 +135,18 @@ runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE + IOobject::NO_WRITE ), linearInterpolate(URec*voidfractionRec) & mesh.Sf() ); - phiRec.write(); + + bool updatePhiRec = false; + if (phiRec.headerOk()) + { + updatePhiRec = true; + phiRec.writeOpt() = IOobject::AUTO_WRITE; + phiRec.write(); + } singlePhaseTransportModel laminarTransport(URec, phiRec); @@ -123,3 +154,40 @@ ( incompressible::turbulenceModel::New(URec, phiRec, laminarTransport) ); + + IOdictionary recDict + ( + IOobject + ( + "recProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + word voidfractionFieldName(recDict.lookupOrDefault("voidfractionFieldName","voidfraction")); + word UFieldName(recDict.lookupOrDefault("UFieldName","U")); + word UsFieldName(recDict.lookupOrDefault("UsFieldName","Us")); + word fluxFieldName(recDict.lookupOrDefault("fluxFieldName","phi")); + + + // place to put weight functions + IOdictionary weightDict + ( + IOobject + ( + "weightDict", + runTime.constant(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ) + ); + if (!weightDict.headerOk()) + { + weightDict.add("weights",scalarList(1,1.0)); + } + scalarList weights(weightDict.lookup("weights")); + Info << "database initial weights: " << weights << endl; diff --git a/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C b/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C index 734e7e8e..2d0c8d6f 100644 --- a/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C +++ b/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C @@ -43,6 +43,7 @@ Rules #include "cfdemCloudRec.H" #include "recBase.H" #include "recModel.H" +#include "recPath.H" #include "cfdemCloud.H" #include "clockModel.H" @@ -59,6 +60,8 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "createFvOptions.H" + #include "readGravitationalAcceleration.H" + cfdemCloudRec particleCloud(mesh); recBase recurrenceBase(mesh); @@ -88,8 +91,9 @@ int main(int argc, char *argv[]) if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 ) { + Info << "updating recurrence fields at time " << runTime.timeName() << "with recTimeIndex = " << recTimeIndex << nl << endl; recurrenceBase.updateRecFields(); - #include "readFields.H" + #include "updateFields.H" recTimeIndex++; } diff --git a/applications/solvers/rcfdemSolverBase/updateFields.H b/applications/solvers/rcfdemSolverBase/updateFields.H new file mode 100644 index 00000000..60a00752 --- /dev/null +++ b/applications/solvers/rcfdemSolverBase/updateFields.H @@ -0,0 +1,29 @@ +scalarList wList(weightDict.lookupOrDefault("weights",scalarList(1,1.0))); + +recurrenceBase.recP().updateIntervalWeights(wList); + +if(recurrenceBase.recM().endOfPath()) +{ + recurrenceBase.extendPath(); +} + +// update fields where necessary +if (updateVoidfractionRec) +{ + recurrenceBase.recM().exportVolScalarField(voidfractionFieldName,voidfractionRec); +} + +if (updateURec) +{ + recurrenceBase.recM().exportVolVectorField(UFieldName,URec); +} + +if (updateUsRec) +{ + recurrenceBase.recM().exportVolVectorField(UsFieldName,UsRec); +} + +if (updatePhiRec) +{ + recurrenceBase.recM().exportSurfaceScalarField(fluxFieldName,phiRec); +} diff --git a/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H b/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H index 8c13a1e4..c79de93c 100644 --- a/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H +++ b/applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H @@ -10,6 +10,7 @@ // main contribution due to gas expansion, not due to transport of kinetic energy // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K) + // assuming constant Cv such that e = Cv * T fvScalarMatrix TEqn = ( fvm::ddt(rhoeps, T) diff --git a/applications/solvers/rcfdemSolverForcedTracers/Make/files b/applications/solvers/rcfdemSolverForcedTracers/Make/files new file mode 100644 index 00000000..c998b354 --- /dev/null +++ b/applications/solvers/rcfdemSolverForcedTracers/Make/files @@ -0,0 +1,3 @@ +rcfdemSolverForcedTracers.C + +EXE=$(CFDEM_APP_DIR)/rcfdemSolverForcedTracers diff --git a/applications/solvers/rcfdemSolverForcedTracers/Make/options b/applications/solvers/rcfdemSolverForcedTracers/Make/options new file mode 100644 index 00000000..9443c2f8 --- /dev/null +++ b/applications/solvers/rcfdemSolverForcedTracers/Make/options @@ -0,0 +1,27 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +EXE_INC = \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ + -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lrecurrence \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools \ + -lfvOptions \ + -l$(CFDEM_LIB_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ +$(CFDEM_ADD_LIBS) diff --git a/applications/solvers/rcfdemSolverForcedTracers/createFields.H b/applications/solvers/rcfdemSolverForcedTracers/createFields.H new file mode 100644 index 00000000..7338925d --- /dev/null +++ b/applications/solvers/rcfdemSolverForcedTracers/createFields.H @@ -0,0 +1,113 @@ + // dummy fields + Info << "\nCreating dummy density field\n" << endl; + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0) + ); + + // particle fields + Info << "\nCreating voidfraction and particle velocity fields\n" << endl; + + volScalarField voidfraction + ( + IOobject + ( + "voidfraction", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // recurrence fields + Info << "\nCreating recurrence fields.\n" << endl; + + volScalarField pRec + ( + IOobject + ( + "pRec", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0) + ); + + volScalarField kRec + ( + IOobject + ( + "kRec", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("k", dimensionSet(0, 2, -2, 0, 0), 0.0) + ); + + volVectorField URec + ( + IOobject + ( + "URec", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//=============================== + + Info << "Calculating face flux field phi\n" << endl; + surfaceScalarField phiRec + ( + IOobject + ( + "phiRec", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(URec*voidfraction) & mesh.Sf() + ); + phiRec.write(); + + singlePhaseTransportModel laminarTransport(URec, phiRec); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(URec, phiRec, laminarTransport) + ); diff --git a/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C b/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C new file mode 100644 index 00000000..8d72b568 --- /dev/null +++ b/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C @@ -0,0 +1,112 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling academic - Open Source CFD-DEM coupling + + Contributing authors: + Thomas Lichtenegger + Copyright (C) 2015- Johannes Kepler University, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling academic. + + CFDEMcoupling academic 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. + + CFDEMcoupling academic 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 CFDEMcoupling academic. If not, see . + +Application + rcfdemSolverForcedTracers + +Description + Moves tracers according to the activated force models on pressure and velocity + fields provided by a recurrence process + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "fvOptions.H" + +#include "recBase.H" +#include "recModel.H" + +#include "cfdemCloud.H" +#include "clockModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFvOptions.H" + + cfdemCloud particleCloud(mesh); + recBase recurrenceBase(mesh); + + const IOdictionary& recProps = mesh.lookupObject("recProperties"); + bool useRecP(recProps.lookupOrDefault("useRecP",false)); + bool useRecK(recProps.lookupOrDefault("useRecK",false)); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl; + + label recTimeIndex = 0; + scalar recTimeStep = recurrenceBase.recM().recTimeStep(); + scalar startTime = runTime.startTime().value(); + + while (runTime.run()) + { + runTime++; + + // do stuff (every lagrangian time step) + particleCloud.clockM().start(1,"Global"); + + Info << "Time = " << runTime.timeName() << nl << endl; + + particleCloud.clockM().start(2,"Coupling"); + + particleCloud.evolve(voidfraction,Us,URec); + + particleCloud.clockM().stop("Coupling"); + + + if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 ) + { + recurrenceBase.updateRecFields(); + #include "updateFields.H" + recTimeIndex++; + } + + particleCloud.clockM().start(27,"Output"); + runTime.write(); + particleCloud.clockM().stop("Output"); + + particleCloud.clockM().stop("Global"); + + Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + } + + Info << "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/rcfdemSolverForcedTracers/updateFields.H b/applications/solvers/rcfdemSolverForcedTracers/updateFields.H new file mode 100644 index 00000000..329710f8 --- /dev/null +++ b/applications/solvers/rcfdemSolverForcedTracers/updateFields.H @@ -0,0 +1,13 @@ +recurrenceBase.recM().exportVolVectorField("U",URec); + +if (useRecP) +{ + recurrenceBase.recM().exportVolScalarField("p",pRec); +} + +if (useRecK) +{ + recurrenceBase.recM().exportVolScalarField("k",kRec); +// in case database contains the velocity variance instead of k, do +// kRec *= 0.5; +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H index d9b32f48..d4011345 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H +++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H @@ -22,19 +22,33 @@ Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); +// For implict T terms in the energy/enthalpy transport equation, use +// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1. +// This formula is valid for ideal gases with e=e(T) and h=h(T). For +// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction +// terms accounting for pressure variations. fvScalarMatrix EEqn ( fvm::div(phi, he) + addSource - Qsource + - QCoeff*T - fvm::Sp(QCoeff/Cpv, he) - // - fvm::laplacian(voidfractionRec*kf/Cpv,he) + + QCoeff/Cpv*he + - fvc::laplacian(voidfractionRec*thCond,T) - fvm::laplacian(voidfractionRec*thCond/Cpv,he) + + fvc::laplacian(voidfractionRec*thCond/Cpv,he) == fvOptions(rho, he) ); + if (transientEEqn) + { + EEqn += fvm::ddt(rho,voidfractionRec,he); + } + + EEqn.relax(); fvOptions.constrain(EEqn); diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H index 454911d9..164799e8 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H +++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H @@ -168,6 +168,8 @@ Info<< "Reading thermophysical properties\n" << endl; linearInterpolate(rho*U*voidfraction) & mesh.Sf() ); + bool transientEEqn(pimple.dict().lookupOrDefault("transientEEqn",false)); + dimensionedScalar rhoMax ( dimensionedScalar::lookupOrDefault diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H new file mode 100644 index 00000000..576fe6a4 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H @@ -0,0 +1,67 @@ +// 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(); + +// For implict T terms in the energy/enthalpy transport equation, use +// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1. +// This formula is valid for ideal gases with e=e(T) and h=h(T). For +// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction +// terms accounting for pressure variations. + + fvScalarMatrix EEqn + ( + fvm::div(phi, he) + + addSource + - Qsource + - QCoeff*T + - fvm::Sp(QCoeff/Cpv, he) + + QCoeff/Cpv*he + - fvc::laplacian(voidfractionRec*thCond,T) + - fvm::laplacian(voidfractionRec*thCond/Cpv,he) + + fvc::laplacian(voidfractionRec*thCond/Cpv,he) + == + fvOptions(rho, he) + ); + + if (transientEEqn) + { + EEqn += fvm::ddt(rho,voidfractionRec,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..44cb1f6e --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H @@ -0,0 +1,35 @@ +// Solve the Momentum equation +particleCloud.otherForces(fOther); + +fvVectorMatrix UEqn +( + fvm::div(phi, U) + + particleCloud.divVoidfractionTau(U, voidfractionRec) + + fvm::Sp(Ksl,U) + - fOther + == + fvOptions(rho, U) +); + +if (stepcounter%nEveryFlow==0) +{ + 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); + } + + + #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..c5cbaca0 --- /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 + ( + mvConvection->fvmDiv(phi, Yi) + - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi) + == + combustion->R(Yi) + + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter() + + fvOptions(rho, Yi) + ); + + YiEqn.relax(); + + fvOptions.constrain(YiEqn); + + YiEqn.solve(mesh.solver("Yi")); + + Yi.relax(); + + 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] + ROOTVSMALL); + forAll(Y,i) + { + if (i!=inertIndex) + { + volScalarField& Yi = Y[i]; + Yi = Yi/(Yt+ROOTVSMALL); + } + } + } + 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..1556801e --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H @@ -0,0 +1,404 @@ +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(thermo.lookupOrDefault("propagateInertSpecie",true)); + + const word inertSpecie(thermo.lookupOrDefault("inertSpecie","none")); + + const scalar inertLowerBound(thermo.lookupOrDefault("inertLowerBound",0.0)); + + const scalar inertUpperBound(thermo.lookupOrDefault("inertUpperBound",1.0)); + + if (!composition.contains(inertSpecie) && inertSpecie != "none") + { + 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() + ); + + bool transientEEqn(pimple.dict().lookupOrDefault("transientEEqn",false)); + + 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) + ); + + 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..c8e4c697 --- /dev/null +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H @@ -0,0 +1,96 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +if (stepcounter%nEveryFlow==0) +{ + +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..e2dc841c --- /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