diff --git a/applications/solvers/cfdemSolverRhoSimple/EEqn.H b/applications/solvers/cfdemSolverRhoSimple/EEqn.H deleted file mode 100644 index a30bbdd0..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/EEqn.H +++ /dev/null @@ -1,57 +0,0 @@ -// 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), voidfraction*U), - p, - "div(phiv,p)" - ) - : fvc::div(phi, K) - ); - - 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; - - fvScalarMatrix EEqn - ( - fvm::div(phi, he) - + addSource - - Qsource - - fvm::Sp(QCoeff/Cpv, he) - - fvm::laplacian(voidfraction*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/cfdemSolverRhoSimple/Make/files b/applications/solvers/cfdemSolverRhoSimple/Make/files deleted file mode 100644 index d9db0744..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -cfdemSolverRhoSimple.C - -EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSimple diff --git a/applications/solvers/cfdemSolverRhoSimple/Make/options b/applications/solvers/cfdemSolverRhoSimple/Make/options deleted file mode 100644 index 0377ece5..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/Make/options +++ /dev/null @@ -1,32 +0,0 @@ -include $(CFDEM_ADD_LIBS_DIR)/additionalLibs - -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$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ - -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ - -EXE_LIBS = \ - -L$(CFDEM_LIB_DIR)\ - -lcompressibleTransportModels \ - -lfluidThermophysicalModels \ - -lspecie \ - -lturbulenceModels \ - -lcompressibleTurbulenceModels \ - -lfiniteVolume \ - -lmeshTools \ - -lsampling \ - -lfvOptions \ - -l$(CFDEM_LIB_COMP_NAME) \ - $(CFDEM_ADD_LIB_PATHS) \ - $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/cfdemSolverRhoSimple/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H deleted file mode 100644 index 8fbd30f2..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/UEqn.H +++ /dev/null @@ -1,30 +0,0 @@ -// Solve the Momentum equation -particleCloud.otherForces(fOther); - -tmp tUEqn -( - fvm::div(phi, U) - + particleCloud.divVoidfractionTau(U, voidfraction) - + 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*Us); -} -else -{ - solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us); -} - -fvOptions.correct(U); - -K = 0.5*magSqr(U); diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C deleted file mode 100644 index 3af7f1f6..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C +++ /dev/null @@ -1,140 +0,0 @@ -/*---------------------------------------------------------------------------*\ -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 - cfdemSolverRhoSimple - -Description - Steady-state solver for turbulent flow of compressible fluids based on - rhoSimpleFoam where functionality for CFD-DEM coupling has been added. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "psiThermo.H" -#include "turbulentFluidThermoModel.H" -#include "bound.H" -#include "simpleControl.H" -#include "fvOptions.H" -#include "localEulerDdtScheme.H" -#include "fvcSmooth.H" - -#include "cfdemCloudEnergy.H" -#include "implicitCouple.H" -#include "clockModel.H" -#include "smoothingModel.H" -#include "forceModel.H" -#include "thermCondModel.H" -#include "energyModel.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" - cfdemCloudEnergy particleCloud(mesh); - #include "checkModelType.H" - - turbulence->validate(); - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (simple.loop()) - { - particleCloud.clockM().start(1,"Global"); - - Info<< "Time = " << runTime.timeName() << nl << endl; - - // do particle stuff - particleCloud.clockM().start(2,"Coupling"); - bool hasEvolved = particleCloud.evolve(voidfraction,Us,U); - - if(hasEvolved) - { - particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces()); - } - - Info << "update Ksl.internalField()" << endl; - Ksl = particleCloud.momCoupleM(0).impMomSource(); - Ksl.correctBoundaryConditions(); - - //Force Checks - vector fTotal(0,0,0); - vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef())); - reduce(fImpTotal, sumOp()); - Info << "TotalForceExp: " << fTotal << endl; - Info << "TotalForceImp: " << fImpTotal << endl; - - #include "solverDebugInfo.H" - particleCloud.clockM().stop("Coupling"); - - particleCloud.clockM().start(26,"Flow"); - - volScalarField rhoeps("rhoeps",rho*voidfraction); - // Pressure-velocity SIMPLE corrector - - #include "UEqn.H" - - - // besides this pEqn, OF offers a "simple consistent"-option - #include "pEqn.H" - rhoeps=rho*voidfraction; - - #include "EEqn.H" - - turbulence->correct(); - - particleCloud.clockM().start(32,"postFlow"); - if(hasEvolved) particleCloud.postFlow(); - particleCloud.clockM().stop("postFlow"); - - 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/cfdemSolverRhoSimple/createFieldRefs.H b/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H deleted file mode 100644 index 5842906a..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H +++ /dev/null @@ -1,2 +0,0 @@ -const volScalarField& T = thermo.T(); -const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/cfdemSolverRhoSimple/createFields.H b/applications/solvers/cfdemSolverRhoSimple/createFields.H deleted file mode 100644 index ab79229b..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/createFields.H +++ /dev/null @@ -1,243 +0,0 @@ -Info<< "Reading thermophysical properties\n" << endl; - - autoPtr pThermo - ( - psiThermo::New(mesh) - ); - psiThermo& thermo = pThermo(); - thermo.validate(args.executable(), "h", "e"); - volScalarField& p = thermo.p(); - - 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 addSource - ( - IOobject - ( - "addSource", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - 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", - simple.dict(), - dimDensity, - GREAT - ) - ); - - dimensionedScalar rhoMin - ( - dimensionedScalar::lookupOrDefault - ( - "rhoMin", - simple.dict(), - dimDensity, - 0 - ) - ); - - Info<< "Creating turbulence model\n" << endl; - autoPtr turbulence - ( - compressible::turbulenceModel::New - ( - rho, - U, - phi, - thermo - ) - ); - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, simple.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)); - - Info<< "\nReading momentum exchange field Ksl\n" << endl; - volScalarField Ksl - ( - IOobject - ( - "Ksl", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0) - ); - - - Info<< "Reading particle velocity field Us\n" << endl; - volVectorField Us - ( - IOobject - ( - "Us", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - -//=============================== diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H deleted file mode 100644 index ed546344..00000000 --- a/applications/solvers/cfdemSolverRhoSimple/pEqn.H +++ /dev/null @@ -1,81 +0,0 @@ -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(voidfraction); -} -volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); - -surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf()); - - -if (simple.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); - - while (simple.correctNonOrthogonal()) - { - // Pressure corrector - fvScalarMatrix pEqn - ( - fvc::div(phi) - - fvm::laplacian(rhorAUf, p) - == - fvOptions(psi, p, rho.name()) - ); - - pEqn.setReference(pRefCell, pRefValue); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi += pEqn.flux(); - } - } -} - -// Explicitly relax pressure for momentum corrector -p.relax(); - -// 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*(voidfraction*fvc::grad(p)-Ksl*Us); -} -else -{ - U = HbyA - rAU*(fvc::grad(p)-Ksl*Us); -} - -U.correctBoundaryConditions(); -fvOptions.correct(U); -K = 0.5*magSqr(U);