mirror of
https://github.com/ParticulateFlow/CFDEMcoupling-PFM.git
synced 2025-12-08 06:37:44 +00:00
Piso with passive temperature equation.
This commit is contained in:
Thomas Lichtenegger
@ -0,0 +1,3 @@
|
||||
cfdemSolverPisoTemp.C
|
||||
|
||||
EXE = $(CFDEM_APP_DIR)/cfdemSolverPisoTemp
|
||||
@ -0,0 +1,19 @@
|
||||
include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
|
||||
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
|
||||
-I$(LIB_SRC)/transportModels \
|
||||
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
|
||||
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
|
||||
|
||||
EXE_LIBS = \
|
||||
-L$(CFDEM_LIB_DIR)\
|
||||
-lincompressibleRASModels \
|
||||
-lincompressibleLESModels \
|
||||
-lincompressibleTransportModels \
|
||||
-lfiniteVolume \
|
||||
-l$(CFDEM_LIB_NAME) \
|
||||
$(CFDEM_ADD_LIB_PATHS) \
|
||||
$(CFDEM_ADD_LIBS)
|
||||
@ -0,0 +1,15 @@
|
||||
// get scalar source from DEM
|
||||
particleCloud.energyContributions(Qsource);
|
||||
Qsource.correctBoundaryConditions();
|
||||
|
||||
// solve temperature transport equation assuming const. density and heat capacity
|
||||
fvScalarMatrix TEqn
|
||||
(
|
||||
fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
|
||||
+ fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
|
||||
- fvm::laplacian(voidfraction*particleCloud.thermCondM().thermDiff(), T)
|
||||
==
|
||||
Qsource/(rho*particleCloud.energyM(0).Cp())
|
||||
);
|
||||
TEqn.relax();
|
||||
TEqn.solve();
|
||||
@ -0,0 +1,201 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
Application
|
||||
cfdemSolverPisoScalar
|
||||
|
||||
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.
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "singlePhaseTransportModel.H"
|
||||
#include "turbulenceModel.H"
|
||||
|
||||
#include "cfdemCloudEnergy.H"
|
||||
#include "implicitCouple.H"
|
||||
#include "smoothingModel.H"
|
||||
#include "forceModel.H"
|
||||
#include "energyModel.H"
|
||||
#include "thermCondModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
#include "setRootCase.H"
|
||||
#include "createTime.H"
|
||||
#include "createMesh.H"
|
||||
#include "createFields.H"
|
||||
#include "initContinuityErrs.H"
|
||||
|
||||
// create cfdemCloud
|
||||
#include "readGravitationalAcceleration.H"
|
||||
cfdemCloudEnergy particleCloud(mesh);
|
||||
#include "checkModelType.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
while (runTime.loop())
|
||||
{
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
#include "readPISOControls.H"
|
||||
#include "CourantNo.H"
|
||||
|
||||
// do particle stuff
|
||||
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();
|
||||
|
||||
|
||||
#include "solverDebugInfo.H"
|
||||
|
||||
#include "TEqn.H"
|
||||
|
||||
if(particleCloud.solveFlow())
|
||||
{
|
||||
// Pressure-velocity PISO corrector
|
||||
{
|
||||
// Momentum predictor
|
||||
fvVectorMatrix UEqn
|
||||
(
|
||||
fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
|
||||
+ fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
|
||||
// + turbulence->divDevReff(U)
|
||||
+ particleCloud.divVoidfractionTau(U, voidfraction)
|
||||
==
|
||||
- fvm::Sp(Ksl/rho,U)
|
||||
);
|
||||
|
||||
UEqn.relax();
|
||||
if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
|
||||
solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
|
||||
else if (momentumPredictor)
|
||||
solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
|
||||
|
||||
// --- PISO loop
|
||||
|
||||
//for (int corr=0; corr<nCorr; corr++)
|
||||
int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
|
||||
|
||||
for (int corr=0; corr<nCorrSoph; corr++)
|
||||
{
|
||||
volScalarField rUA = 1.0/UEqn.A();
|
||||
|
||||
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
|
||||
volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
|
||||
surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
|
||||
|
||||
U = rUA*UEqn.H();
|
||||
|
||||
#ifdef version23
|
||||
phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
|
||||
+ rUAfvoidfraction*fvc::ddtCorr(U, phi);
|
||||
#else
|
||||
phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
|
||||
+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
|
||||
#endif
|
||||
surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
|
||||
surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
|
||||
|
||||
if (modelType=="A")
|
||||
rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
// Pressure corrector
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
|
||||
);
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
if
|
||||
(
|
||||
corr == nCorr-1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
{
|
||||
pEqn.solve(mesh.solver("pFinal"));
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve();
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phiGes -= pEqn.flux();
|
||||
phi = phiGes;
|
||||
}
|
||||
|
||||
} // end non-orthogonal corrector loop
|
||||
|
||||
#include "continuityErrorPhiPU.H"
|
||||
|
||||
if (modelType=="B" || modelType=="Bfull")
|
||||
U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
|
||||
else
|
||||
U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
|
||||
|
||||
U.correctBoundaryConditions();
|
||||
|
||||
} // end piso loop
|
||||
}
|
||||
|
||||
turbulence->correct();
|
||||
}// end solveFlow
|
||||
else
|
||||
{
|
||||
Info << "skipping flow solution." << endl;
|
||||
}
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,156 @@
|
||||
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
|
||||
);
|
||||
|
||||
//========================
|
||||
// drag law modelling
|
||||
//========================
|
||||
|
||||
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(0, 0, -1, 0, 0), 1.0)
|
||||
);
|
||||
|
||||
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::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("0", dimensionSet(1, -3, 0, 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
|
||||
);
|
||||
|
||||
//========================
|
||||
// scalar field modelling
|
||||
//========================
|
||||
|
||||
Info<< "\nCreating temperature field\n" << endl;
|
||||
volScalarField T
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"T",
|
||||
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::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*voidfraction) & mesh.Sf()
|
||||
);
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
|
||||
|
||||
|
||||
singlePhaseTransportModel laminarTransport(U, phi);
|
||||
|
||||
autoPtr<incompressible::turbulenceModel> turbulence
|
||||
(
|
||||
incompressible::turbulenceModel::New(U, phi, laminarTransport)
|
||||
);
|
||||
Reference in New Issue
Block a user