Merge pull request #123 from tmjnijssen/develop

TUe/SMR developments
This commit is contained in:
Daniel
2021-11-05 09:07:27 +01:00
committed by GitHub
181 changed files with 837983 additions and 235 deletions

View File

@ -55,14 +55,21 @@ Foam::multiphaseMixture::calcNu() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
// 1/nu
tmp<volScalarField> tnuInv = iter()/iter().nu();
volScalarField& nuInv = tnuInv.ref();
// nu
tmp<volScalarField> tnu = iter()*iter().nu(); tmp<volScalarField> tnu = iter()*iter().nu();
volScalarField& nu = tnu.ref(); volScalarField& nu = tnu.ref();
for (++iter; iter != phases_.end(); ++iter) for (++iter; iter != phases_.end(); ++iter)
{ {
nu += iter()*iter().nu(); nuInv += iter()/iter().nu();
} }
nu = 1/nuInv;
return tnu; return tnu;
} }

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
set -x
wclean libso multiphaseMixture
wclean
#------------------------------------------------------------------------------

View File

@ -0,0 +1,12 @@
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Parse arguments for library compilation
targetType=libso
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
set -x
wmake $targetType multiphaseMixture
wmake
#------------------------------------------------------------------------------

View File

@ -0,0 +1,22 @@
// get mixture properties
Cs = mixture.Cs();
diffusionCorrection = mixture.diffusionCorrection();
Deff = particleCloud.diffCoeffM().diffCoeff();
// get scalar source from DEM
particleCloud.massContributions(Sm);
particleCloud.massCoefficients(Smi);
fvScalarMatrix CEqn
(
fvm::ddt(voidfraction,C)
+ fvm::div(phi,C)
- fvm::laplacian(Deff*voidfraction,C)
+ fvm::div(fvc::interpolate(Deff*voidfraction)*diffusionCorrection*mesh.magSf(), C)
==
Sm + fvm::Sp(Smi,C)
);
CEqn.relax();
fvOptions.constrain(CEqn);
CEqn.solve();

View File

@ -0,0 +1,22 @@
// get mixture properties
Cp = mixture.Cp();
kf = mixture.kf();
// get scalar source from DEM
particleCloud.energyContributions(Qsource);
particleCloud.energyCoefficients(QCoeff);
fvScalarMatrix EEqn
(
rho*Cp*(fvm::ddt(voidfraction,T)
+ fvm::div(phi,T))
- fvm::laplacian(thCond*voidfraction,T)
==
Qsource + fvm::Sp(QCoeff,T)
);
EEqn.relax();
fvOptions.constrain(EEqn);
EEqn.solve();

View File

@ -0,0 +1,3 @@
cfdemSolverMultiphaseScalar.C
EXE = $(CFDEM_APP_DIR)/cfdemSolverMultiphaseScalar

View File

@ -0,0 +1,35 @@
FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
EXE_INC = \
$(PFLAGS) \
-I$(CFDEM_OFVERSION_DIR) \
-ImultiphaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-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 \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\
-lcfdemMultiphaseInterFoamScalar \
-linterfaceProperties \
-lincompressibleTransportModels \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-lsampling \
-l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS)

View File

@ -0,0 +1,61 @@
const surfaceScalarField& rhoPhi(mixture.rhoPhi());
volScalarField muEff = rho*(turbulence->nu() + turbulence->nut());
if (modelType == "A")
muEff *= voidfraction;
fvVectorMatrix UEqn
(
fvm::ddt(rhoEps, U) - fvm::Sp(fvc::ddt(rhoEps),U)
+ fvm::div(rhoPhi, U) - fvm::Sp(fvc::div(rhoPhi),U)
//+ particleCloud.divVoidfractionTau(U, voidfraction)
- fvm::laplacian(muEff, U) - fvc::div(muEff*dev2(fvc::grad(U)().T()))
==
fvOptions(rho, U)
- fvm::Sp(Ksl,U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
{
solve
(
UEqn
==
fvc::reconstruct
(
(- ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh)) * mesh.magSf()
)
+
fvc::reconstruct
(
mixture.surfaceTensionForce() * mesh.magSf()
) * voidfraction
+ Ksl*Us
);
fvOptions.correct(U);
}
else if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
) * voidfraction
+ Ksl*Us
);
fvOptions.correct(U);
}

View File

@ -0,0 +1,17 @@
// Additional solver-specific checks
// Useful if one wants to e.g. initialize floating particles using the Archimedes model
if (particleCloud.couplingProperties().found("unrestrictedForceModelSelection"))
{
Warning << "Using unrestrictedForceModelSelection, results may be incorrect!" << endl;
} else
{
#include "checkModelType.H"
}
word modelType = particleCloud.modelType();
if(!particleCloud.couplingProperties().found("useDDTvoidfraction"))
{
Warning << "Suppressing ddt(voidfraction) is not recommended with this solver as it may generate incorrect results!" << endl;
}

View File

@ -0,0 +1,21 @@
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi
(
mixture.nearInterface()().primitiveField()
*fvc::surfaceSum(mag(phi))().primitiveField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanAlphaCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

View File

@ -0,0 +1,164 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
Application
cfdemSolverMultiphaseScalar
Description
CFD-DEM solver for n incompressible fluids which captures the interfaces and
includes surface-tension and contact-angle effects for each phase. It is based
on the OpenFOAM(R)-4.x solver multiphaseInterFoam but extended to incorporate
DEM functionalities from the open-source DEM code LIGGGHTS.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "multiphaseMixture.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "cfdemCloudEnergy.H"
#include "implicitCouple.H"
#include "clockModel.H"
#include "smoothingModel.H"
#include "forceModel.H"
#include "thermCondModel.H"
#include "diffCoeffModel.H"
#include "energyModel.H"
#include "massTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#if OPENFOAM_VERSION_MAJOR >= 6
FatalError << "cfdemSolverMultiphase requires OpenFOAM 4.x or 5.x to work properly" << exit(FatalError);
#endif
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "correctPhi.H"
#include "CourantNo.H"
turbulence->validate();
// create cfdemCloud
cfdemCloudEnergy particleCloud(mesh);
#include "additionalChecks.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
particleCloud.clockM().start(1,"Global");
Info<< "Time = " << runTime.timeName() << nl << endl;
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.internalField()*(Us.internalField()-U.internalField())).value();
reduce(fImpTotal, sumOp<vector>());
Info << "TotalForceExp: " << fTotal << endl;
Info << "TotalForceImp: " << fImpTotal << endl;
#include "solverDebugInfo.H"
particleCloud.clockM().stop("Coupling");
particleCloud.clockM().start(26,"Flow");
if(particleCloud.solveFlow())
{
mixture.solve();
rho = mixture.rho();
rhoEps = rho * voidfraction;
#include "EEqn.H"
#include "CEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
}
else
{
Info << "skipping flow solution." << endl;
}
particleCloud.clockM().start(31,"postFlow");
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;
}
// ************************************************************************* //

View File

@ -0,0 +1,11 @@
CorrectPhi
(
U,
phi,
p_rgh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
geometricZeroField(),
pimple
);
#include "continuityErrs.H"

View File

@ -0,0 +1,342 @@
//===============================
// particle interaction 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(1, -3, -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
);
voidfraction.oldTime();
Info<< "Reading particle velocity field Us\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
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()
);
multiphaseMixture mixture(U, phi, voidfraction);
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mixture.rho()
);
rho.oldTime();
//========================
// scalar field modelling
//========================
Info<< "Reading/creating thermal fields\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
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)
);
volScalarField QCoeff
(
IOobject
(
"Qsource",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
);
volScalarField Cp
(
IOobject
(
"Cp",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mixture.Cp()
);
volScalarField kf
(
IOobject
(
"kf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mixture.kf()
);
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<< "Reading/creating concentration fields\n" << endl;
volScalarField C
(
IOobject
(
"C",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField Sm
(
IOobject
(
"Sm",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
);
volScalarField Smi
(
IOobject
(
"Smi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(0,0,-1,0,0,0,0), 0.0)
);
volScalarField D
(
IOobject
(
"D",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mixture.D()
);
volScalarField Deff
(
IOobject
(
"Deff",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
);
volScalarField Cs
(
IOobject
(
"Cs",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mixture.Cs()
);
surfaceScalarField diffusionCorrection
(
IOobject
(
"diffusionCorrection",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mixture.diffusionCorrection()
);
//========================
volScalarField rhoEps ("rhoEps", rho * voidfraction);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}
mesh.setFluxRequired(p_rgh.name());

View File

@ -0,0 +1,5 @@
phase/phase.C
alphaContactAngle/alphaContactAngleFvPatchScalarField.C
multiphaseMixture.C
LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoamScalar

View File

@ -0,0 +1,18 @@
FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
EXE_INC = \
$(PFLAGS) \
-IalphaContactAngle \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-Wno-deprecated-copy
LIB_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools

View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "alphaContactAngleFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
(
Istream& is
)
:
theta0_(readScalar(is)),
uTheta_(readScalar(is)),
thetaA_(readScalar(is)),
thetaR_(readScalar(is))
{}
Istream& operator>>
(
Istream& is,
alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
)
{
is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
return is;
}
Ostream& operator<<
(
Ostream& os,
const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
)
{
os << tp.theta0_ << token::SPACE
<< tp.uTheta_ << token::SPACE
<< tp.thetaA_ << token::SPACE
<< tp.thetaR_;
return os;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
zeroGradientFvPatchScalarField(p, iF)
{}
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const alphaContactAngleFvPatchScalarField& gcpsf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
thetaProps_(gcpsf.thetaProps_)
{}
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
zeroGradientFvPatchScalarField(p, iF),
thetaProps_(dict.lookup("thetaProperties"))
{
evaluate();
}
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const alphaContactAngleFvPatchScalarField& gcpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
zeroGradientFvPatchScalarField(gcpsf, iF),
thetaProps_(gcpsf.thetaProps_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
os.writeKeyword("thetaProperties")
<< thetaProps_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
alphaContactAngleFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,215 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::alphaContactAngleFvPatchScalarField
Description
Contact-angle boundary condition for multi-phase interface-capturing
simulations. Used in conjuction with multiphaseMixture.
SourceFiles
alphaContactAngleFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef alphaContactAngleFvPatchScalarField_H
#define alphaContactAngleFvPatchScalarField_H
#include "zeroGradientFvPatchFields.H"
#include "multiphaseMixture.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class alphaContactAngleFvPatch Declaration
\*---------------------------------------------------------------------------*/
class alphaContactAngleFvPatchScalarField
:
public zeroGradientFvPatchScalarField
{
public:
class interfaceThetaProps
{
//- Equilibrium contact angle
scalar theta0_;
//- Dynamic contact angle velocity scale
scalar uTheta_;
//- Limiting advancing contact angle
scalar thetaA_;
//- Limiting receeding contact angle
scalar thetaR_;
public:
// Constructors
interfaceThetaProps()
{}
interfaceThetaProps(Istream&);
// Member functions
//- Return the equilibrium contact angle theta0
scalar theta0(bool matched=true) const
{
if (matched) return theta0_;
else return 180.0 - theta0_;
}
//- Return the dynamic contact angle velocity scale
scalar uTheta() const
{
return uTheta_;
}
//- Return the limiting advancing contact angle
scalar thetaA(bool matched=true) const
{
if (matched) return thetaA_;
else return 180.0 - thetaA_;
}
//- Return the limiting receeding contact angle
scalar thetaR(bool matched=true) const
{
if (matched) return thetaR_;
else return 180.0 - thetaR_;
}
// IO functions
friend Istream& operator>>(Istream&, interfaceThetaProps&);
friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
};
typedef HashTable
<
interfaceThetaProps,
multiphaseMixture::interfacePair,
multiphaseMixture::interfacePair::hash
> thetaPropsTable;
private:
// Private data
thetaPropsTable thetaProps_;
public:
//- Runtime type information
TypeName("alphaContactAngle");
// Constructors
//- Construct from patch and internal field
alphaContactAngleFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
alphaContactAngleFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given alphaContactAngleFvPatchScalarField
// onto a new patch
alphaContactAngleFvPatchScalarField
(
const alphaContactAngleFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField>
(
new alphaContactAngleFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
alphaContactAngleFvPatchScalarField
(
const alphaContactAngleFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new alphaContactAngleFvPatchScalarField(*this, iF)
);
}
// Member functions
//- Return the contact angle properties
const thetaPropsTable& thetaProps() const
{
return thetaProps_;
}
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,929 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "multiphaseMixture.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "Time.H"
#include "subCycle.H"
#include "MULES.H"
#include "surfaceInterpolate.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
#include "fvcDiv.H"
#include "fvcFlux.H"
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
const Foam::scalar Foam::multiphaseMixture::convertToRad =
Foam::constant::mathematical::pi/180.0;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::multiphaseMixture::calcAlphas()
{
scalar level = 0.0;
alphas_ == 0.0;
forAllIter(PtrDictionary<phase>, phases_, iter)
{
alphas_ += level*iter();
level += 1.0;
}
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::calcNu() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
// 1/nu
tmp<volScalarField> tnuInv = iter()/iter().nu();
volScalarField& nuInv = tnuInv.ref();
// nu
tmp<volScalarField> tnu = iter()*iter().nu();
volScalarField& nu = tnu.ref();
for (++iter; iter != phases_.end(); ++iter)
{
nuInv += iter()/iter().nu();
}
nu = 1/nuInv;
return tnu;
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::calcStf() const
{
tmp<surfaceScalarField> tstf
(
new surfaceScalarField
(
IOobject
(
"stf",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar
(
"stf",
dimensionSet(1, -2, -2, 0, 0),
0.0
)
)
);
surfaceScalarField& stf = tstf.ref();
forAllConstIter(PtrDictionary<phase>, phases_, iter1)
{
const phase& alpha1 = iter1();
PtrDictionary<phase>::const_iterator iter2 = iter1;
++iter2;
for (; iter2 != phases_.end(); ++iter2)
{
const phase& alpha2 = iter2();
sigmaTable::const_iterator sigma =
sigmas_.find(interfacePair(alpha1, alpha2));
if (sigma == sigmas_.end())
{
FatalErrorInFunction
<< "Cannot find interface " << interfacePair(alpha1, alpha2)
<< " in list of sigma values"
<< exit(FatalError);
}
stf += dimensionedScalar("sigma", dimSigma_, sigma())
*fvc::interpolate(K(alpha1, alpha2))*
(
fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
- fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
);
}
}
return tstf;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::multiphaseMixture::multiphaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& voidfraction
)
:
IOdictionary
(
IOobject
(
"transportProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
phases_(lookup("phases"), phase::iNew(U, phi)),
mesh_(U.mesh()),
U_(U),
phi_(phi),
voidfraction_(voidfraction),
rhoPhi_
(
IOobject
(
"rhoPhi",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
),
surfaceTensionForce_
(
IOobject
(
"surfaceTensionForce",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar("surfaceTensionForce", dimensionSet(1, -2, -2, 0, 0), 0.0)
),
alphas_
(
IOobject
(
"alphas",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar("alphas", dimless, 0.0)
),
nu_
(
IOobject
(
"nu",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
calcNu()
),
sigmas_(lookup("sigmas")),
dimSigma_(1, 0, -2, 0, 0),
deltaN_
(
"deltaN",
1e-8/pow(average(mesh_.V()), 1.0/3.0)
)
{
calcAlphas();
alphas_.write();
surfaceTensionForce_ = calcStf();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::rho() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
tmp<volScalarField> trho = iter()*iter().rho();
volScalarField& rho = trho.ref();
for (++iter; iter != phases_.end(); ++iter)
{
rho += iter()*iter().rho();
}
return trho;
}
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::rho(const label patchi) const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
scalarField& rho = trho.ref();
for (++iter; iter != phases_.end(); ++iter)
{
rho += iter().boundaryField()[patchi]*iter().rho().value();
}
return trho;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::mu() const
{
Info << "In multiphasemixture mu()" << endl;
return rho()*nu();
// PtrDictionary<phase>::const_iterator iter = phases_.begin();
// tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
// volScalarField& mu = tmu.ref();
// for (++iter; iter != phases_.end(); ++iter)
// {
// mu += iter()*iter().rho()*iter().nu();
// }
// return tmu;
}
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::mu(const label patchi) const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
tmp<scalarField> tmu =
iter().boundaryField()[patchi]
*iter().rho().value()
*iter().nu(patchi);
scalarField& mu = tmu.ref();
for (++iter; iter != phases_.end(); ++iter)
{
mu +=
iter().boundaryField()[patchi]
*iter().rho().value()
*iter().nu(patchi);
}
return tmu;
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::muf() const
{
return nuf()*fvc::interpolate(rho());
// PtrDictionary<phase>::const_iterator iter = phases_.begin();
// tmp<surfaceScalarField> tmuf =
// fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
// surfaceScalarField& muf = tmuf.ref();
// for (++iter; iter != phases_.end(); ++iter)
// {
// muf +=
// fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
// }
// return tmuf;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nu() const
{
return nu_;
}
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::nu(const label patchi) const
{
//return nu_.boundaryField()[patchi];
PtrDictionary<phase>::const_iterator iter = phases_.begin();
tmp<scalarField> tnu =
iter().boundaryField()[patchi]
*iter().nu(patchi);
scalarField& nu = tnu.ref();
for (++iter; iter != phases_.end(); ++iter)
{
nu +=
iter().boundaryField()[patchi]
*iter().nu(patchi);
}
return tnu;
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::nuf() const
{
//return muf()/fvc::interpolate(rho());
PtrDictionary<phase>::const_iterator iter = phases_.begin();
tmp<surfaceScalarField> tnuf =
fvc::interpolate(iter())*fvc::interpolate(iter().nu());
surfaceScalarField& nuf = tnuf.ref();
for (++iter; iter != phases_.end(); ++iter)
{
nuf +=
fvc::interpolate(iter())*fvc::interpolate(iter().nu());
}
return tnuf;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cp() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
// rho*Cp
tmp<volScalarField> trhoCp = iter()*iter().Cp()*iter().rho();
volScalarField& rhoCp = trhoCp.ref();
// Cp
tmp<volScalarField> tCp = iter()*iter().Cp();
volScalarField& Cp = tCp.ref();
for (++iter; iter != phases_.end(); ++iter)
{
rhoCp += iter()*iter().Cp()*iter().rho();
}
Cp = rhoCp/rho();
return tCp;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::kf() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
// rho*Cp/kf
tmp<volScalarField> trhoCpkf = iter()*iter().rho()*iter().Cp()/iter().kf();
volScalarField& rhoCpkf = trhoCpkf.ref();
// kf
tmp<volScalarField> tkf = iter()*iter().kf();
volScalarField& kf = tkf.ref();
for (++iter; iter != phases_.end(); ++iter)
{
rhoCpkf += iter()*iter().rho()*iter().Cp()/iter().kf();
}
kf = rho()*Cp()/rhoCpkf;
return tkf;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::D() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
// 1/D
tmp<volScalarField> tDInv = iter()/iter().D();
volScalarField& DInv = tDInv.ref();
// D
tmp<volScalarField> tD = iter()*iter().D();
volScalarField& D = tD.ref();
for (++iter; iter != phases_.end(); ++iter)
{
DInv += iter()/iter().D();
}
D = 1/DInv;
return tD;
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cs() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
// Cs
tmp<volScalarField> tCs = iter()*iter().Cs();
volScalarField& Cs = tCs.ref();
for (++iter; iter != phases_.end(); ++iter)
{
Cs += iter()*iter().Cs();
}
return tCs;
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::diffusionCorrection() const
{
surfaceScalarField numerator
(
IOobject
(
"numerator",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimless/dimLength, 0.0)
);
surfaceScalarField denominator
(
IOobject
(
"denominator",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimless, 0.0)
);
PtrDictionary<phase>::const_iterator iter = phases_.begin();
const phase& alpha1 = iter();
for (++iter; iter != phases_.end(); ++iter)
{
const phase& alpha2 = iter();
scalar He = alpha1.Cs().value() / (alpha2.Cs().value() + SMALL);
numerator += (1/He - 1) * fvc::snGrad(alpha2);
denominator += fvc::interpolate(alpha2) * (1/He - 1);
}
tmp<surfaceScalarField> correction = numerator / (denominator + 1 + SMALL);
/*
PtrDictionary<phase>::const_iterator iter = phases_.begin();
const phase& alphaL = iter();
++iter;
const phase& alphaG = iter();
scalar He = alphaG.Cs().value() / (alphaL.Cs().value() + SMALL);
surfaceScalarField gradAlphaL = fvc::snGrad(alphaL);
surfaceScalarField surfAlphaL = fvc::interpolate(alphaL);
tmp<surfaceScalarField> correction = (1-He)/(surfAlphaL + He*(1-surfAlphaL) + 10*SMALL) * gradAlphaL;
*/
return correction;
}
void Foam::multiphaseMixture::solve()
{
correct();
const Time& runTime = mesh_.time();
volScalarField& alpha = phases_.first();
const dictionary& alphaControls = mesh_.solverDict("alpha");
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh_
),
mesh_,
dimensionedScalar("0", rhoPhi_.dimensions(), 0)
);
dimensionedScalar totalDeltaT = runTime.deltaT();
for
(
subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
FatalError << "Sub-cycling of the alpha equation not yet implemented!!" << abort(FatalError);
solveAlphas(cAlpha);
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
}
rhoPhi_ = rhoPhiSum;
}
else
{
solveAlphas(cAlpha);
}
// Update the mixture kinematic viscosity
nu_ = calcNu();
surfaceTensionForce_ = calcStf();
}
void Foam::multiphaseMixture::correct()
{
forAllIter(PtrDictionary<phase>, phases_, iter)
{
iter().correct();
}
}
Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
(
const volScalarField& alpha1,
const volScalarField& alpha2
) const
{
/*
// Cell gradient of alpha
volVectorField gradAlpha =
alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
// Interpolated face-gradient of alpha
surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
*/
surfaceVectorField gradAlphaf
(
fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
- fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
);
// Face unit interface normal
return gradAlphaf/(mag(gradAlphaf) + deltaN_);
}
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
(
const volScalarField& alpha1,
const volScalarField& alpha2
) const
{
// Face unit interface normal flux
return nHatfv(alpha1, alpha2) & mesh_.Sf();
}
// Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.
// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.
void Foam::multiphaseMixture::correctContactAngle
(
const phase& alpha1,
const phase& alpha2,
surfaceVectorField::Boundary& nHatb
) const
{
const volScalarField::Boundary& gbf
= alpha1.boundaryField();
const fvBoundaryMesh& boundary = mesh_.boundary();
forAll(boundary, patchi)
{
if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
{
const alphaContactAngleFvPatchScalarField& acap =
refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
vectorField& nHatPatch = nHatb[patchi];
vectorField AfHatPatch
(
mesh_.Sf().boundaryField()[patchi]
/mesh_.magSf().boundaryField()[patchi]
);
alphaContactAngleFvPatchScalarField::thetaPropsTable::
const_iterator tp =
acap.thetaProps().find(interfacePair(alpha1, alpha2));
if (tp == acap.thetaProps().end())
{
FatalErrorInFunction
<< "Cannot find interface " << interfacePair(alpha1, alpha2)
<< "\n in table of theta properties for patch "
<< acap.patch().name()
<< exit(FatalError);
}
bool matched = (tp.key().first() == alpha1.name());
scalar theta0 = convertToRad*tp().theta0(matched);
scalarField theta(boundary[patchi].size(), theta0);
scalar uTheta = tp().uTheta();
// Calculate the dynamic contact angle if required
if (uTheta > SMALL)
{
scalar thetaA = convertToRad*tp().thetaA(matched);
scalar thetaR = convertToRad*tp().thetaR(matched);
// Calculated the component of the velocity parallel to the wall
vectorField Uwall
(
U_.boundaryField()[patchi].patchInternalField()
- U_.boundaryField()[patchi]
);
Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
// Find the direction of the interface parallel to the wall
vectorField nWall
(
nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
);
// Normalise nWall
nWall /= (mag(nWall) + SMALL);
// Calculate Uwall resolved normal to the interface parallel to
// the interface
scalarField uwall(nWall & Uwall);
theta += (thetaA - thetaR)*tanh(uwall/uTheta);
}
// Reset nHatPatch to correspond to the contact angle
scalarField a12(nHatPatch & AfHatPatch);
scalarField b1(cos(theta));
scalarField b2(nHatPatch.size());
forAll(b2, facei)
{
b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
scalarField det(1.0 - a12*a12);
scalarField a((b1 - a12*b2)/det);
scalarField b((b2 - a12*b1)/det);
nHatPatch = a*AfHatPatch + b*nHatPatch;
nHatPatch /= (mag(nHatPatch) + deltaN_.value());
}
}
}
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
(
const phase& alpha1,
const phase& alpha2
) const
{
tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
// Simple expression for curvature
return -fvc::div(tnHatfv & mesh_.Sf());
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nearInterface() const
{
tmp<volScalarField> tnearInt
(
new volScalarField
(
IOobject
(
"nearInterface",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("nearInterface", dimless, 0.0)
)
);
forAllConstIter(PtrDictionary<phase>, phases_, iter)
{
tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
}
return tnearInt;
}
void Foam::multiphaseMixture::solveAlphas
(
const scalar cAlpha
)
{
static label nSolves=-1;
nSolves++;
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phic(mag(phi_/mesh_.magSf()));
phic = min(cAlpha*phic, max(phic));
PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
int phasei = 0;
forAllIter(PtrDictionary<phase>, phases_, iter)
{
phase& alpha = iter();
alphaPhiCorrs.set
(
phasei,
new surfaceScalarField
(
"phi" + alpha.name() + "Corr",
fvc::flux
(
phi_,
alpha,
alphaScheme
)
)
);
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
forAllIter(PtrDictionary<phase>, phases_, iter2)
{
phase& alpha2 = iter2();
if (&alpha2 == &alpha) continue;
surfaceScalarField phir(phic*nHatf(alpha, alpha2));
alphaPhiCorr += fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha,
alpharScheme
);
}
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
voidfraction_,
alpha,
phi_,
alphaPhiCorr,
zeroField(),
zeroField(),
#if OPENFOAM_VERSION_MAJOR < 6
1,
0,
#else
oneField(),
zeroField(),
#endif
true
);
phasei++;
}
MULES::limitSum(alphaPhiCorrs);
rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
volScalarField sumAlpha
(
IOobject
(
"sumAlpha",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("sumAlpha", dimless, 0)
);
phasei = 0;
forAllIter(PtrDictionary<phase>, phases_, iter)
{
phase& alpha = iter();
surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
MULES::explicitSolve
(
voidfraction_,
alpha,
alphaPhi,
zeroField(),
zeroField()
);
rhoPhi_ += alphaPhi*alpha.rho();
Info<< alpha.name() << " volume fraction, min, max = "
<< alpha.weightedAverage(mesh_.V()).value()
<< ' ' << min(alpha).value()
<< ' ' << max(alpha).value()
<< endl;
sumAlpha += alpha;
phasei++;
}
Info<< "Phase-sum volume fraction, min, max = "
<< sumAlpha.weightedAverage(mesh_.V()).value()
<< ' ' << min(sumAlpha).value()
<< ' ' << max(sumAlpha).value()
<< endl;
calcAlphas();
}
bool Foam::multiphaseMixture::read()
{
if (transportModel::read())
{
bool readOK = true;
PtrList<entry> phaseData(lookup("phases"));
label phasei = 0;
forAllIter(PtrDictionary<phase>, phases_, iter)
{
readOK &= iter().read(phaseData[phasei++].dict());
}
lookup("sigmas") >> sigmas_;
return readOK;
}
else
{
return false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,299 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
Class
multiphaseMixture
Description
This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
which is an incompressible multi-phase mixture with built in solution
for the phase fractions with interface compression for interface-capturing.
It has been extended to include the void fraction in the volume fraction
transport equations.
Derived from transportModel so that it can be unsed in conjunction with
the incompressible turbulence models.
Surface tension and contact-angle is handled for the interface
between each phase-pair.
SourceFiles
multiphaseMixture.C
\*---------------------------------------------------------------------------*/
#ifndef multiphaseMixture_H
#define multiphaseMixture_H
#include "incompressible/transportModel/transportModel.H"
#include "IOdictionary.H"
#include "phase.H"
#include "PtrDictionary.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class multiphaseMixture Declaration
\*---------------------------------------------------------------------------*/
class multiphaseMixture
:
public IOdictionary,
public transportModel
{
public:
class interfacePair
:
public Pair<word>
{
public:
class hash
:
public Hash<interfacePair>
{
public:
hash()
{}
label operator()(const interfacePair& key) const
{
return word::hash()(key.first()) + word::hash()(key.second());
}
};
// Constructors
interfacePair()
{}
interfacePair(const word& alpha1Name, const word& alpha2Name)
:
Pair<word>(alpha1Name, alpha2Name)
{}
interfacePair(const phase& alpha1, const phase& alpha2)
:
Pair<word>(alpha1.name(), alpha2.name())
{}
// Friend Operators
friend bool operator==
(
const interfacePair& a,
const interfacePair& b
)
{
return
(
((a.first() == b.first()) && (a.second() == b.second()))
|| ((a.first() == b.second()) && (a.second() == b.first()))
);
}
friend bool operator!=
(
const interfacePair& a,
const interfacePair& b
)
{
return (!(a == b));
}
};
private:
// Private data
//- Dictionary of phases
PtrDictionary<phase> phases_;
const fvMesh& mesh_;
const volVectorField& U_;
const surfaceScalarField& phi_;
const volScalarField& voidfraction_;
surfaceScalarField rhoPhi_;
surfaceScalarField surfaceTensionForce_;
volScalarField alphas_;
volScalarField nu_;
typedef HashTable<scalar, interfacePair, interfacePair::hash>
sigmaTable;
sigmaTable sigmas_;
dimensionSet dimSigma_;
//- Stabilisation for normalisation of the interface normal
const dimensionedScalar deltaN_;
//- Conversion factor for degrees into radians
static const scalar convertToRad;
// Private member functions
void calcAlphas();
tmp<volScalarField> calcNu() const;
void solveAlphas(const scalar cAlpha);
tmp<surfaceVectorField> nHatfv
(
const volScalarField& alpha1,
const volScalarField& alpha2
) const;
tmp<surfaceScalarField> nHatf
(
const volScalarField& alpha1,
const volScalarField& alpha2
) const;
void correctContactAngle
(
const phase& alpha1,
const phase& alpha2,
surfaceVectorField::Boundary& nHatb
) const;
tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;
tmp<surfaceScalarField> calcStf() const;
public:
// Constructors
//- Construct from components
multiphaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& voidfraction
);
//- Destructor
virtual ~multiphaseMixture()
{}
// Member Functions
//- Return the phases
const PtrDictionary<phase>& phases() const
{
return phases_;
}
//- Return the velocity
const volVectorField& U() const
{
return U_;
}
//- Return the volumetric flux
const surfaceScalarField& phi() const
{
return phi_;
}
const surfaceScalarField& rhoPhi() const
{
return rhoPhi_;
}
//- Return the mixture density
tmp<volScalarField> rho() const;
//- Return the mixture density for patch
tmp<scalarField> rho(const label patchi) const;
//- Return the dynamic laminar viscosity
tmp<volScalarField> mu() const;
//- Return the dynamic laminar viscosity for patch
tmp<scalarField> mu(const label patchi) const;
//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> muf() const;
//- Return the kinematic laminar viscosity
tmp<volScalarField> nu() const;
//- Return the laminar viscosity for patch
tmp<scalarField> nu(const label patchi) const;
//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> nuf() const;
//- Return the heat capacity
tmp<volScalarField> Cp() const;
//- Return the thermal conductivity
tmp<volScalarField> kf() const;
//- Return the diffusion coefficient
tmp<volScalarField> D() const;
//- Return the solubility
tmp<volScalarField> Cs() const;
//- Return the diffusion correction term
tmp<surfaceScalarField> diffusionCorrection() const;
tmp<surfaceScalarField> surfaceTensionForce() const
{
return surfaceTensionForce_;
}
//- Indicator of the proximity of the interface
// Field values are 1 near and 0 away for the interface.
tmp<volScalarField> nearInterface() const;
//- Solve for the mixture phase-fractions
void solve();
//- Correct the mixture properties
void correct();
//- Read base transportProperties dictionary
bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "phase.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::phase::phase
(
const word& phaseName,
const dictionary& phaseDict,
const volVectorField& U,
const surfaceScalarField& phi
)
:
volScalarField
(
IOobject
(
IOobject::groupName("alpha", phaseName),
U.mesh().time().timeName(),
U.mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
U.mesh()
),
name_(phaseName),
phaseDict_(phaseDict),
nuModel_
(
viscosityModel::New
(
IOobject::groupName("nu", phaseName),
phaseDict_,
U,
phi
)
),
rho_("rho", dimDensity, phaseDict_),
Cp_("Cp", (dimSpecificHeatCapacity), phaseDict_),
kf_("kf", (dimPower/dimLength/dimTemperature), phaseDict_),
D_("D", dimViscosity, phaseDict_),
Cs_("Cs", dimDensity, phaseDict_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::phase> Foam::phase::clone() const
{
NotImplemented;
return autoPtr<phase>(NULL);
}
void Foam::phase::correct()
{
nuModel_->correct();
}
bool Foam::phase::read(const dictionary& phaseDict)
{
phaseDict_ = phaseDict;
phaseDict_.lookup("Cp") >> Cp_;
phaseDict_.lookup("kf") >> kf_;
phaseDict_.lookup("D") >> D_;
phaseDict_.lookup("Cs") >> Cs_;
if (nuModel_->read(phaseDict_))
{
phaseDict_.lookup("rho") >> rho_;
return true;
}
else
{
return false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,190 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::phase
Description
Single incompressible phase derived from the phase-fraction.
Used as part of the multiPhaseMixture for interface-capturing multi-phase
simulations.
SourceFiles
phase.C
\*---------------------------------------------------------------------------*/
#ifndef phase_H
#define phase_H
#include "volFields.H"
#include "dictionaryEntry.H"
#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phase Declaration
\*---------------------------------------------------------------------------*/
class phase
:
public volScalarField
{
// Private data
word name_;
dictionary phaseDict_;
autoPtr<viscosityModel> nuModel_;
dimensionedScalar rho_;
dimensionedScalar Cp_;
dimensionedScalar kf_;
dimensionedScalar D_;
dimensionedScalar Cs_;
public:
// Constructors
//- Construct from components
phase
(
const word& name,
const dictionary& phaseDict,
const volVectorField& U,
const surfaceScalarField& phi
);
//- Return clone
autoPtr<phase> clone() const;
//- Return a pointer to a new phase created on freestore
// from Istream
class iNew
{
const volVectorField& U_;
const surfaceScalarField& phi_;
public:
iNew
(
const volVectorField& U,
const surfaceScalarField& phi
)
:
U_(U),
phi_(phi)
{}
autoPtr<phase> operator()(Istream& is) const
{
dictionaryEntry ent(dictionary::null, is);
return autoPtr<phase>(new phase(ent.keyword(), ent, U_, phi_));
}
};
// Member Functions
const word& name() const
{
return name_;
}
const word& keyword() const
{
return name();
}
//- Return const-access to phase1 viscosityModel
const viscosityModel& nuModel() const
{
return nuModel_();
}
//- Return the kinematic laminar viscosity
tmp<volScalarField> nu() const
{
return nuModel_->nu();
}
//- Return the laminar viscosity for patch
tmp<scalarField> nu(const label patchi) const
{
return nuModel_->nu(patchi);
}
//- Return const-access to phase1 density
const dimensionedScalar& rho() const
{
return rho_;
}
//- Return const-access to phase1 heat capacity
const dimensionedScalar& Cp() const
{
return Cp_;
}
//- Return const-access to phase1 thermal conductivity
const dimensionedScalar& kf() const
{
return kf_;
}
//- Return const-access to phase1 diffusion coefficient
const dimensionedScalar& D() const
{
return D_;
}
//- Return const-access to phase1 solubility
const dimensionedScalar& Cs() const
{
return Cs_;
}
//- Correct the phase properties
void correct();
//-Inherit read from volScalarField
using volScalarField::read;
//- Read base transportProperties dictionary
bool read(const dictionary& phaseDict);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,73 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUepsf("rAUepsf", fvc::interpolate(rAU*voidfraction));
surfaceScalarField rAUepsSqf("rAUepsSqf", fvc::interpolate(rAU*voidfraction*voidfraction));
volVectorField Ueps("Ueps", U * voidfraction);
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA*voidfraction)
+ fvc::interpolate(voidfraction*rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
if (modelType == "A")
rAUepsf = rAUepsSqf;
surfaceScalarField phig (-ghf*fvc::snGrad(rho)*rAUepsf*mesh.magSf());
surfaceScalarField phiSt (mixture.surfaceTensionForce()*rAUepsSqf*mesh.magSf());
surfaceScalarField phiS (fvc::flux(voidfraction*Us*Ksl*rAU));
phiHbyA += phig + phiSt + phiS;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, Ueps, phiHbyA, rAUepsf);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUepsf, p_rgh) == particleCloud.ddtVoidfraction() + fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();
p_rgh.relax();
if (modelType == "A")
U = HbyA + voidfraction*rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
else
U = HbyA + rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -73,7 +73,8 @@ compressible, reacting flows.
heatTransferGranConduction, heatTransferGranConduction,
heatTransferGunn, heatTransferGunn,
heatTransferRanzMarshall, heatTransferRanzMarshall,
reactionHeat :tb(c=2,ea=c) reactionHeat,
wallHeatTransferYagi :tb(c=2,ea=c)
6.7 Force models :h4 6.7 Force models :h4
@ -113,6 +114,7 @@ particleDeformation,
potentialRelaxation, potentialRelaxation,
"surfaceTensionForce"_forceModel_surfaceTensionForce.html, "surfaceTensionForce"_forceModel_surfaceTensionForce.html,
"virtualMassForce"_forceModel_virtualMassForce.html, "virtualMassForce"_forceModel_virtualMassForce.html,
"ParmarBassetForce"_forceModel_ParmarBassetForce.html,
"viscForce"_forceModel_viscForce.html, "viscForce"_forceModel_viscForce.html,
"volWeightedAverage"_forceModel_volWeightedAverage.html :tb(c=2,ea=c) "volWeightedAverage"_forceModel_volWeightedAverage.html :tb(c=2,ea=c)
@ -203,7 +205,8 @@ smoothing the exchange fields.
"constDiffSmoothing"_smoothingModel_constDiffSmoothing.html, "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html,
"off"_smoothingModel_noSmoothing.html, "off"_smoothingModel_noSmoothing.html,
"temporalSmoothing"_smoothingModel_temporalSmoothing.html :tb(c=2,ea=c) "temporalSmoothing"_smoothingModel_temporalSmoothing.html,
"constDiffAndTemporalSmoothing"_smoothingModel_constDiffAndTemporalSmoothing.html :tb(c=2,ea=c)
6.16 Thermal conductivity models :h4 6.16 Thermal conductivity models :h4
@ -229,3 +232,18 @@ accounting for the volume of the particles in the CFD domain.
off, off,
trilinear :tb(c=2,ea=c) trilinear :tb(c=2,ea=c)
6.18 Mass transfer models :h4
The {massTransferModels} keyword specifies a list of mass transfer models used evaluating
species transfer between particles and fluids.
massTransferGunn :tb(c=2,ea=c)
6.19 Diffusion coefficient models :h4
The {diffCoeffModel} keyword entry specifies the model for the diffusion
coefficient of dissolved spieces in the fluid phase in the presence of particles.
SyamlalDiffCoeff,
off :tb(c=2,ea=c)

View File

@ -11,6 +11,7 @@ This section lists all CFDEMcoupling solvers alphabetically.
"cfdemSolverIB"_cfdemSolverIB.html, "cfdemSolverIB"_cfdemSolverIB.html,
"cfdemSolverMultiphase"_cfdemSolverMultiphase.html, "cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
"cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html,
"cfdemSolverPiso"_cfdemSolverPiso.html, "cfdemSolverPiso"_cfdemSolverPiso.html,
"cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html, "cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html,
"cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html, "cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html,

View File

@ -0,0 +1,68 @@
<!-- HTML_ONLY -->
<HEAD>
<META CHARSET="utf-8">
</HEAD>
<!-- END_HTML_ONLY -->
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
cfdemSolverMultiphaseScalar command :h3
[Description:]
<!-- HTML_ONLY -->
"cfdemSolverMultiphaseScalar" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam&reg;(*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles. Additionally, it enable evaluation of temperature and dissolved species concentration fields.
<!-- END_HTML_ONLY -->
<!-- RST
"cfdemSolverMultiphaseScalar" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam\ |reg|\ (*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles. Additionally, it enable evaluation of temperature and dissolved species concentration fields.
.. |reg| unicode:: U+000AE .. REGISTERED SIGN
END_RST -->
For more details, see "Vångö et al. (2018)"_#Vångö2018 and "Nijssen et al. (2021)"_#Nijssen2021.
:line
:link(Vångö2018)
[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018):
"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds",
Applied Mathematical Modelling
:link(Nijssen2021)
[(Nijssen2021)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2021):
"article title here",
journal name here
:line
<!-- HTML_ONLY -->
NOTE:
(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com, and owner of the
OPENFOAM&reg; and OpenCFD&reg; trade marks.
OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com.
<!-- END_HTML_ONLY -->
<!-- RST
.. note::
(*) This offering is not approved or endorsed by OpenCFD Limited, producer
and distributor of the OpenFOAM software via www.openfoam.com, and owner of
the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com.
.. |reg| unicode:: U+000AE .. REGISTERED SIGN
END_RST -->

View File

@ -0,0 +1,83 @@
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
energyModel wallHeatTransferYagi command :h3
[Syntax:]
Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
dictionary.
energyModels
(
wallHeatTransferYagi
);
wallHeatTransferYagiProps
\{
QWallFluidName "QWallFluid";
QWallFluidCoeffNAme "QWallFluidCoeff";
wallTempName "wallTemp";
tempFieldName "T";
voidfractionFieldName "voidfraction";
voidfractionMax scalar1;
maxSource scalar2;
velFieldName "U";
densityFieldName "rho";
implicit switch1;
verbose switch2;
\} :pre
{QWallFluid} = name of the wall-fluid heat transfer rate field :ulb,l
{QwallFluidCoeff} = name of the wall-fluid heat transfer coefficient field :l
{wallTemp} = name of the wall temperature field :l
{T} = name of the fluid temperature field :l
{voidfraction} = name of the finite volume void fraction field :l
{scalar1} = maximum void fraction for a cell to be considered packed :l
{scalar2} = (optional, default 1e30) maximum allowed source term :l
{U} = name of the finite volume fluid velocity field :l
{rho} = name of the fluid density field :l
{switch1} = (optional, default true) activate to treat the fluid temperature implicitly :l
{switch2} = (optional, default false) activate to write additional fields, and write cell-based data to the terminal :l
:ule
[Examples:]
energyModels
(
wallHeatTransferYagi
);
wallHeatTransferYagiProps
\{
QWallFluidName "QWallFluid";
QWallFluidCoeffNAme "QWallFluidCoeff";
wallTempName "wallTemp";
tempFieldName "T";
voidfractionFieldName "voidfraction";
voidfractionMax 0.5;
velFieldName "U";
densityFieldName "rho";
implicit true;
\} :pre
[Description:]
The energy model performs calculation of wall-to-bed heat transfer for a
packed bed, based on the correlation of "Yagi and Wakao (1959)"_#Yagi1959.
:link(Yagi1959)
[(Yagi1959)] S. Yagi, N. Wakao. (1959):
"Heat and Mass Transfer from Wall to Fluid in Packed Beds",
AIChE Journal
[Restrictions:]
IMPORTANT NOTE: Model not validated!
[Related commands:]
energyModel

View File

@ -19,19 +19,23 @@ forceModels
MeiLiftProps MeiLiftProps
\{ \{
velFieldName "U"; velFieldName "U";
useSecondOrderTerms; useShearInduced switch1;
treatForceExplicit switch1; useSpinInduced switch2;
verbose switch2; combineShearSpin switch3;
interpolation switch3; treatForceExplicit switch4;
scalarViscosity switch4; verbose switch5;
interpolation switch6;
scalarViscosity switch7;
\} :pre \} :pre
{U} = name of the finite volume fluid velocity field :ulb,l {U} = name of the finite volume fluid velocity field :ulb,l
{useSecondOrderTerms} = (optional, default false) switch to activate second order terms in the lift force model :l {useShearInduced} = (optional, default true) switch to activate shear-induced (Saffman) lift :l
{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l {useSpinInduced} = (optional, default false) switch to activate spin-induced (Magnus) lift :l
{switch2} = (optional, default false) switch to activate the report of per-particle quantities to the screen :l {combineShearSpin} = (optional, default false) switch to use the correlation for combined shear- and spin-induced lift by "Loth and Dorgan (2009)"_#LothDorgan2009 :l
{switch3} = (optional, default false) switch to activate tri-linear interpolation of the flow quantities at the particle position :l
{switch4} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l {switch4} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{switch5} = (optional, default false) switch to activate the report of per-particle quantities to the screen :l
{switch6} = (optional, default false) switch to activate tri-linear interpolation of the flow quantities at the particle position :l
{switch7} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
:ule :ule
[Examples:] [Examples:]
@ -43,7 +47,6 @@ forceModels
MeiLiftProps MeiLiftProps
\{ \{
velFieldName "U"; velFieldName "U";
useSecondOrderTerms;
interpolation true; interpolation true;
verbose true; verbose true;
\} :pre \} :pre
@ -52,9 +55,12 @@ MeiLiftProps
The force model performs the calculation of forces (e.g. fluid-particle The force model performs the calculation of forces (e.g. fluid-particle
interaction forces) acting on each DEM particle. The {MeiLift} model calculates interaction forces) acting on each DEM particle. The {MeiLift} model calculates
the lift force for each particle based on Loth and Dorgan (2009). In case the the lift force for each particle based on "Loth and Dorgan (2009)"_#LothDorgan2009.
keyword "useSecondOrderTerms" is not specified, this lift force model uses the
expression of McLaughlin (1991, J. Fluid Mech. 224:261-274). :link(LothDorgan2009)
[(LothDorgan2009)] E. Loth, A.J. Dorgan. (2009):
"An equation of motion for particles of finite Reynolds number and size",
Environmental Fluid Mechanics
[Restrictions:] [Restrictions:]

View File

@ -0,0 +1,88 @@
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
forceModel ParmarBassetForce command :h3
[Syntax:]
Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
dictionary.
forceModels
(
ParmarBassetForce
);
ParmarBassetForceProps
\{
velFieldName "U";
UsFieldName "Us";
nIntegral scalar1;
discretisationOrder scalar2;
treatForceExplicit switch1;
interpolation switch2;
smoothingModel "smoothingModel";
\} :pre
{U} = name of the finite volume fluid velocity field :ulb,l
{Us} = name of the finite volume cell averaged particle velocity field :l
{scalar1} = number of timesteps considered in the near history :l
{scalar2} = (1 or 2) discretisation order of the far history differential equations :l
{switch1} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{smoothingModel} = name of smoothing model for the dU/dt field :l
:ule
[Examples:]
forceModels
(
ParmarBassetForce
);
ParmarBassetForceProps
\{
velFieldName "U";
USFieldName "Us"
nIntegral 10;
discretisationOrder 1;
smoothingModel constDiffAndTemporalSmoothing;
constDiffAndTemporalSmoothingProps
\{
lowerLimit 1e-8;
upperLimit 1e8;
smoothingLength 0.1;
smoothingStrength 0.001;
correctBoundary true;
\}
\} :pre
[Description:]
The force model performs the calculation of forces (e.g. fluid-particle
interaction forces) acting on each DEM particle. The {ParmarBassetForce} model
calculates the Basset history force for each particle, based on the method by
"Parmar et al. (2018)"_#Parmar2018.
For more detail, see "Nijssen et al. (2020)"_#Nijssen2020.
:link(Parmar2018)
[(Parmar2018)] M. Parmar, S. Annamalai, S. Balachandar, A. Prosperetti. (2018):
"Differential Formulation of the Viscous History Force on a Particle for Efficient and Accurate Computation",
Journal of Fluid Mechanics
:link(Nijssen2020)
[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
International Journal of Multiphase Flow
[Restrictions:]
IMPORTANT NOTE: Model not validated!
[Related commands:]
"forceModel"_forceModel.html

View File

@ -19,19 +19,29 @@ forceModels
virtualMassForceProps virtualMassForceProps
\{ \{
velFieldName "U"; velFieldName "U";
voidfractionFieldName "voidfraction";
UsFieldName "Us";
phiFieldName "phi"; phiFieldName "phi";
splitUrelCalculation switch1; splitUrelCalculation switch1;
useUs switch2;
useFelderhof switch3;
Cadd scalar1; Cadd scalar1;
treatForceExplicit switch2; treatForceExplicit switch4;
interpolation switch3; interpolation switch5;
smoothingModel "smoothingModel";
\} :pre \} :pre
{U} = name of the finite volume fluid velocity field :ulb,l {U} = name of the finite volume fluid velocity field :ulb,l
{voidfraction} = name of the finite volume void fraction field :l
{Us} = name of the finite volume cell averaged particle velocity field :l
{phi} = name of the finite volume flux field :l {phi} = name of the finite volume flux field :l
{switch1} = (optional, default false) indicator to split calculation of Urel between CFDEM and LIGGGHTS :l {switch1} = (optional, default false) indicator to split calculation of Urel between CFDEM and LIGGGHTS :l
{switch2} = (optional, default false) indicator to use the cell-averaged particle velocity Us for calculation of Urel :l
{switch3} = (optional, default false) indicator to use the correlation by Felderhof for the virtual mass coefficient :l
{scalar1} = (optional, default 0.5) virtual mass coefficient :l {scalar1} = (optional, default 0.5) virtual mass coefficient :l
{switch2} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l {switch4} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{switch3} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l {switch5} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{smoothingModel} = name of smoothing model for the DU/Dt field :l
:ule :ule
[Examples:] [Examples:]
@ -43,7 +53,10 @@ forceModels
virtualMassForceProps virtualMassForceProps
\{ \{
velFieldName "U"; velFieldName "U";
voidfractionFieldName "voidfraction"
USFieldName "Us"
phiFieldName "phi"; phiFieldName "phi";
smoothingModel off;
\} :pre \} :pre
[Description:] [Description:]
@ -52,6 +65,15 @@ The force model performs the calculation of forces (e.g. fluid-particle
interaction forces) acting on each DEM particle. The {virtualMassForce} model interaction forces) acting on each DEM particle. The {virtualMassForce} model
calculates the virtual mass force for each particle. calculates the virtual mass force for each particle.
For more detail, see "Nijssen et al. (2020)"_#Nijssen2020.
:link(Nijssen2020)
[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
International Journal of Multiphase Flow
[Restrictions:]
IMPORTANT NOTE: Model not validated! IMPORTANT NOTE: Model not validated!
[Related commands:] [Related commands:]

View File

@ -0,0 +1,77 @@
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
smoothingModel constDiffAndTemporalSmoothing command :h3
[Syntax:]
Defined in dictionary depending on the application.
smoothingModel constDiffAndTemporalSmoothing;
constDiffAndTemporalSmoothingProps
\{
lowerLimit number1;
upperLimit number2;
smoothingLength lengthScale;
smoothingLengthReferenceField lengthScaleRefField;
smoothingStrength smoothingStrength;
correctBoundary switch1;
verbose;
\} :pre
{number1} = scalar fields will be bound to this lower value :ulb,l
{number2} = scalar fields will be bound to this upper value :l
{lengthScale} = length scale over which the exchange fields will be smoothed out :l
{lengthScaleRefField} = length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
{smoothingStrength} = control parameter gamma for the smoothing, lower value yields stronger smoothing (gamma = 1 results in an equal contribution from the un-smoothed and smoothed fields) :l
{correctBoundary} = (optional, default false) activate to use purely temporal smoothing on the boundary field, avoids interpolation errors near the domain boundary :l
{verbose} = (optional, default false) flag for debugging output :l
:ule
[Examples:]
constDiffAndTemporalSmoothingProps
\{
lowerLimit 0.1;
upperLimit 1e10;
smoothingLength 1500e-6;
smoothingLengthReferenceField 9000e-6;
referenceField "p";
gamma 1.0;
\} :pre
[Description:]
The {constDiffAndTemporalSmoothing} model is a smoothing model that combines the
spacial smoothing of "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html and
the temporal smoothing of "temporalSmoothing"_smoothingModel_temporalSmoothing.html for
the relaxation of a desired quantity. This model can be used to filter out
high-frequency fluctuations (e.g. numerical noise) controlled via the smoothing length
and the temporal smoothing strength parameter gamma.
For more details, see "Vångö et al. (2018)"_#Vångö2018 and "Nijssen et al. (2020)"_#Nijssen2020.
:link(Vångö2018)
[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018):
"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds",
Applied Mathematical Modelling
:link(Nijssen2020)
[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
International Journal of Multiphase Flow
[Restrictions:]
This model is tested in a limited number of flow situations.
[Related commands:]
"smoothingModel"_smoothingModel.html
"constDiffSmoothing"_smoothingModel_constDiffSmoothing.html
"temporalSmoothing"_smoothingModel_temporalSmoothing.html

View File

@ -221,7 +221,7 @@ alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh
alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh' alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh'
#- shortcut to visualize the clock model data #- shortcut to visualize the clock model data
alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py' alias vizClock='python2 $CFDEM_UT_DIR/vizClock/matPlot.py'
#- recursive touch of current directory #- recursive touch of current directory
alias touchRec='find ./* -exec touch {} \;' alias touchRec='find ./* -exec touch {} \;'

View File

@ -3,3 +3,4 @@ lagrangian/cfdemParticleComp/dir
recurrence/dir recurrence/dir
finiteVolume/dir finiteVolume/dir
../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir ../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir
../applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/dir

View File

@ -12,6 +12,7 @@ cfdemSolverIB/dir
cfdemSolverPisoScalar/dir cfdemSolverPisoScalar/dir
cfdemSolverRhoPimpleChem/dir cfdemSolverRhoPimpleChem/dir
cfdemSolverMultiphase/dir cfdemSolverMultiphase/dir
cfdemSolverMultiphaseScalar/dir
rcfdemSolverRhoSteadyPimpleChem/dir rcfdemSolverRhoSteadyPimpleChem/dir
rctfSpeciesTransport/dir rctfSpeciesTransport/dir
cfdemSolverPisoFreeStreaming/dir cfdemSolverPisoFreeStreaming/dir

View File

@ -48,7 +48,8 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
densityFieldName_("rho"), densityFieldName_("rho"),
tubeLength_(-1.), tubeLength_(-1.),
tubeDiameter_(-1.), tubeDiameter_(-1.),
zeta_(-1) zeta_(-1),
p0_(-1)
{ {
Info << "error, wrong constructor1!" << endl; Info << "error, wrong constructor1!" << endl;
} }
@ -71,7 +72,8 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
densityFieldName_("rho"), densityFieldName_("rho"),
tubeLength_(ptf.tubeLength_), tubeLength_(ptf.tubeLength_),
tubeDiameter_(ptf.tubeDiameter_), tubeDiameter_(ptf.tubeDiameter_),
zeta_(ptf.zeta_) zeta_(ptf.zeta_),
p0_(ptf.p0_)
{ {
const scalar t = this->db().time().timeOutputValue(); const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t)); fvPatchField<Type>::operator==(uniformValue_->value(t));
@ -89,13 +91,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
: :
fixedValueFvPatchField<Type>(p, iF), fixedValueFvPatchField<Type>(p, iF),
uniformValue_(Function1<Type>::New("uniformValue", dict)), uniformValue_(Function1<Type>::New("uniformValue", dict)),
pName_("p"), //JOKER pName_("p_rgh"), //JOKER
phiName_("phi"), //JOKER phiName_("phi"), //JOKER
velocityFieldName_("U"), velocityFieldName_("U"),
densityFieldName_("rho"), densityFieldName_("rho"),
tubeLength_(readScalar(dict.lookup("tubeLength"))), tubeLength_(readScalar(dict.lookup("tubeLength"))),
tubeDiameter_(readScalar(dict.lookup("tubeDiameter"))), tubeDiameter_(readScalar(dict.lookup("tubeDiameter"))),
zeta_(readScalar(dict.lookup("zeta"))) zeta_(readScalar(dict.lookup("zeta"))),
p0_(readScalar(dict.lookup("p0")))
{ {
const scalar t = this->db().time().timeOutputValue(); const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t)); fvPatchField<Type>::operator==(uniformValue_->value(t));
@ -110,13 +113,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
: :
fixedValueFvPatchField<Type>(ptf), fixedValueFvPatchField<Type>(ptf),
uniformValue_(ptf.uniformValue_().clone().ptr()), uniformValue_(ptf.uniformValue_().clone().ptr()),
pName_("p"), //JOKER pName_("p_rgh"), //JOKER
phiName_("phi"), //JOKER phiName_("phi"), //JOKER
velocityFieldName_("U"), velocityFieldName_("U"),
densityFieldName_("rho"), densityFieldName_("rho"),
tubeLength_(ptf.tubeLength_), tubeLength_(ptf.tubeLength_),
tubeDiameter_(ptf.tubeDiameter_), tubeDiameter_(ptf.tubeDiameter_),
zeta_(ptf.zeta_) zeta_(ptf.zeta_),
p0_(ptf.p0_)
{ {
const scalar t = this->db().time().timeOutputValue(); const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t)); fvPatchField<Type>::operator==(uniformValue_->value(t));
@ -132,13 +136,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
: :
fixedValueFvPatchField<Type>(ptf, iF), fixedValueFvPatchField<Type>(ptf, iF),
uniformValue_(ptf.uniformValue_().clone().ptr()), uniformValue_(ptf.uniformValue_().clone().ptr()),
pName_("p"), //JOKER pName_("p_rgh"), //JOKER
phiName_("phi"), //JOKER phiName_("phi"), //JOKER
velocityFieldName_("U"), velocityFieldName_("U"),
densityFieldName_("rho"), densityFieldName_("rho"),
tubeLength_(ptf.tubeLength_), tubeLength_(ptf.tubeLength_),
tubeDiameter_(ptf.tubeDiameter_), tubeDiameter_(ptf.tubeDiameter_),
zeta_(ptf.zeta_) zeta_(ptf.zeta_),
p0_(ptf.p0_)
{ {
const scalar t = this->db().time().timeOutputValue(); const scalar t = this->db().time().timeOutputValue();
fvPatchField<Type>::operator==(uniformValue_->value(t)); fvPatchField<Type>::operator==(uniformValue_->value(t));
@ -185,7 +190,7 @@ void uniformFixedValueTubeFvPatchField<Type>::updateCoeffs()
// calc pressure drop // calc pressure drop
// fvPatchField<Type>::operator==(0.5*zeta*density*mag(velocity)*mag(velocity)*uniformValue_->value(t)); // fvPatchField<Type>::operator==(0.5*zeta*density*mag(velocity)*mag(velocity)*uniformValue_->value(t));
//dP = zeta * l/d * rho u^2 / 2 //dP = zeta * l/d * rho u^2 / 2
fvPatchField<Type>::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()-pressure)*uRelFact)*uniformValue_->value(t)); fvPatchField<Type>::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()+p0_-pressure)*uRelFact)*uniformValue_->value(t));
//fvPatchField<Type>::operator==((100.*uniformValue_->value(t))); //fvPatchField<Type>::operator==((100.*uniformValue_->value(t)));
fixedValueFvPatchField<Type>::updateCoeffs(); fixedValueFvPatchField<Type>::updateCoeffs();
} }
@ -200,6 +205,7 @@ void uniformFixedValueTubeFvPatchField<Type>::write(Ostream& os) const
os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl; os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl;
os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl; os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl;
os.writeKeyword("zeta") << zeta_ << token::END_STATEMENT << nl; os.writeKeyword("zeta") << zeta_ << token::END_STATEMENT << nl;
os.writeKeyword("p0") << p0_ << token::END_STATEMENT << nl;
} }

View File

@ -74,6 +74,8 @@ class uniformFixedValueTubeFvPatchField
scalar zeta_; scalar zeta_;
scalar p0_;
public: public:
//- Runtime type information //- Runtime type information

View File

@ -2,10 +2,12 @@ cfdemCloud = cfdemCloud
cfdTools = cfdTools cfdTools = cfdTools
recModels = subModels/recModel recModels = subModels/recModel
energyModels = subModels/energyModel energyModels = subModels/energyModel
massTransferModels = subModels/massTransferModel
forceModels = subModels/forceModel forceModels = subModels/forceModel
forceSubModels = subModels/forceModel/forceSubModels forceSubModels = subModels/forceModel/forceSubModels
forceModelsMS = subModels/forceModelMS forceModelsMS = subModels/forceModelMS
thermCondModels = subModels/thermCondModel thermCondModels = subModels/thermCondModel
diffCoeffModels = subModels/diffCoeffModel
chemistryModels = subModels/chemistryModel chemistryModels = subModels/chemistryModel
IOModels = subModels/IOModel IOModels = subModels/IOModel
voidFractionModels = subModels/voidFractionModel voidFractionModels = subModels/voidFractionModel
@ -44,14 +46,24 @@ $(energyModels)/energyModel/newEnergyModel.C
$(energyModels)/heatTransferGunn/heatTransferGunn.C $(energyModels)/heatTransferGunn/heatTransferGunn.C
$(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C $(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C
$(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C $(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C
$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C
$(energyModels)/reactionHeat/reactionHeat.C $(energyModels)/reactionHeat/reactionHeat.C
$(massTransferModels)/massTransferModel/massTransferModel.C
$(massTransferModels)/massTransferModel/newMassTransferModel.C
$(massTransferModels)/massTransferGunn/massTransferGunn.C
$(thermCondModels)/thermCondModel/thermCondModel.C $(thermCondModels)/thermCondModel/thermCondModel.C
$(thermCondModels)/thermCondModel/newThermCondModel.C $(thermCondModels)/thermCondModel/newThermCondModel.C
$(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C $(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
$(thermCondModels)/noTherm/noThermCond.C $(thermCondModels)/noTherm/noThermCond.C
$(diffCoeffModels)/diffCoeffModel/diffCoeffModel.C
$(diffCoeffModels)/diffCoeffModel/newDiffCoeffModel.C
$(diffCoeffModels)/SyamlalDiffCoeff/SyamlalDiffCoeff.C
$(diffCoeffModels)/noDiff/noDiffCoeff.C
$(forceModels)/forceModel/forceModel.C $(forceModels)/forceModel/forceModel.C
$(forceModels)/forceModel/newForceModel.C $(forceModels)/forceModel/newForceModel.C
$(forceModels)/noDrag/noDrag.C $(forceModels)/noDrag/noDrag.C
@ -68,6 +80,7 @@ $(forceModels)/KochHillDrag/KochHillDrag.C
$(forceModels)/KochHillRWDrag/KochHillRWDrag.C $(forceModels)/KochHillRWDrag/KochHillRWDrag.C
$(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
$(forceModels)/virtualMassForce/virtualMassForce.C $(forceModels)/virtualMassForce/virtualMassForce.C
$(forceModels)/ParmarBassetForce/ParmarBassetForce.C
$(forceModels)/gradPForce/gradPForce.C $(forceModels)/gradPForce/gradPForce.C
$(forceModels)/viscForce/viscForce.C $(forceModels)/viscForce/viscForce.C
$(forceModels)/MeiLift/MeiLift.C $(forceModels)/MeiLift/MeiLift.C
@ -191,5 +204,6 @@ $(smoothingModels)/smoothingModel/newSmoothingModel.C
$(smoothingModels)/noSmoothing/noSmoothing.C $(smoothingModels)/noSmoothing/noSmoothing.C
$(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
$(smoothingModels)/temporalSmoothing/temporalSmoothing.C $(smoothingModels)/temporalSmoothing/temporalSmoothing.C
$(smoothingModels)/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C
LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME) LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

View File

@ -82,6 +82,7 @@ cfdemCloud::cfdemCloud
solveFlow_(couplingProperties_.lookupOrDefault<bool>("solveFlow", true)), solveFlow_(couplingProperties_.lookupOrDefault<bool>("solveFlow", true)),
verbose_(couplingProperties_.found("verbose")), verbose_(couplingProperties_.found("verbose")),
ignore_(couplingProperties_.found("ignore")), ignore_(couplingProperties_.found("ignore")),
multiphase_(couplingProperties_.lookupOrDefault<bool>("multiphase",false)),
allowCFDsubTimestep_(true), allowCFDsubTimestep_(true),
modelCheck_(couplingProperties_.lookupOrDefault<bool>("modelCheck",true)), modelCheck_(couplingProperties_.lookupOrDefault<bool>("modelCheck",true)),
limitDEMForces_(couplingProperties_.found("limitDEMForces")), limitDEMForces_(couplingProperties_.found("limitDEMForces")),
@ -89,6 +90,7 @@ cfdemCloud::cfdemCloud
getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)), getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)), getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)), getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
getParticleAngVels_(couplingProperties_.lookupOrDefault<bool>("getParticleAngVels",false)),
streamingMode_(couplingProperties_.lookupOrDefault<bool>("streamingMode",false)), streamingMode_(couplingProperties_.lookupOrDefault<bool>("streamingMode",false)),
streamingFluc_(couplingProperties_.lookupOrDefault<bool>("streamingFluc",false)), streamingFluc_(couplingProperties_.lookupOrDefault<bool>("streamingFluc",false)),
maxDEMForce_(0.), maxDEMForce_(0.),
@ -108,6 +110,7 @@ cfdemCloud::cfdemCloud
particleDensities_(NULL), particleDensities_(NULL),
particleEffVolFactors_(NULL), particleEffVolFactors_(NULL),
particleTypes_(NULL), particleTypes_(NULL),
particleAngVels_(NULL),
particleWeights_(NULL), particleWeights_(NULL),
particleVolumes_(NULL), particleVolumes_(NULL),
particleV_(NULL), particleV_(NULL),
@ -406,6 +409,7 @@ cfdemCloud::~cfdemCloud()
if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1); if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1); if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1); if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
if(getParticleAngVels_) dataExchangeM().destroy(particleAngVels_,1);
for for
( (
@ -437,6 +441,7 @@ void cfdemCloud::getDEMdata()
if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_); if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_);
if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_); if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_);
if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_); if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_);
if(getParticleAngVels_) dataExchangeM().getData("omega","vector-atom",particleAngVels_);
} }
void cfdemCloud::giveDEMdata() void cfdemCloud::giveDEMdata()
@ -493,7 +498,12 @@ void cfdemCloud::setForces()
resetArray(expForces_,numberOfParticles(),3); resetArray(expForces_,numberOfParticles(),3);
resetArray(DEMForces_,numberOfParticles(),3); resetArray(DEMForces_,numberOfParticles(),3);
resetArray(Cds_,numberOfParticles(),1); resetArray(Cds_,numberOfParticles(),1);
for (int i=0;i<cfdemCloud::nrForceModels();i++) cfdemCloud::forceM(i).setForce(); for (int i=0;i<cfdemCloud::nrForceModels();i++)
{
clockM().start(30+i,forceModels_[i]);
cfdemCloud::forceM(i).setForce();
clockM().stop(forceModels_[i]);
}
if (limitDEMForces_) if (limitDEMForces_)
{ {
@ -811,6 +821,7 @@ bool cfdemCloud::reAllocArrays()
if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1); if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1); if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1); if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
if(getParticleAngVels_) dataExchangeM().allocateArray(particleAngVels_,0.,3);
for for
( (

View File

@ -93,6 +93,8 @@ protected:
const bool ignore_; const bool ignore_;
const bool multiphase_;
bool allowCFDsubTimestep_; bool allowCFDsubTimestep_;
const bool modelCheck_; const bool modelCheck_;
@ -107,6 +109,8 @@ protected:
const bool getParticleTypes_; const bool getParticleTypes_;
const bool getParticleAngVels_;
const bool streamingMode_; const bool streamingMode_;
const bool streamingFluc_; const bool streamingFluc_;
@ -145,6 +149,8 @@ protected:
int **particleTypes_; // per particle type int **particleTypes_; // per particle type
double **particleAngVels_; // particle angular velocities
double **particleWeights_; // particle weights double **particleWeights_; // particle weights
double **particleVolumes_; // particle volumes with subcell information double **particleVolumes_; // particle volumes with subcell information
@ -321,6 +327,8 @@ public:
inline bool verbose() const; inline bool verbose() const;
inline bool multiphase() const;
inline const IOdictionary& couplingProperties() const; inline const IOdictionary& couplingProperties() const;
inline double ** positions() const; inline double ** positions() const;
@ -382,6 +390,10 @@ public:
virtual inline int ** particleTypes() const; virtual inline int ** particleTypes() const;
virtual inline label particleType(label index) const; virtual inline label particleType(label index) const;
inline bool getParticleAngVels() const;
virtual inline double ** particleAngVels() const;
virtual vector particleAngVel(label index) const;
//access to the particle's rotation and torque data //access to the particle's rotation and torque data
virtual inline double ** DEMTorques() const {return NULL;} virtual inline double ** DEMTorques() const {return NULL;}
virtual inline double ** omegaArray() const {return NULL;} virtual inline double ** omegaArray() const {return NULL;}

View File

@ -128,6 +128,11 @@ inline bool cfdemCloud::verbose() const
return verbose_; return verbose_;
} }
inline bool cfdemCloud::multiphase() const
{
return multiphase_;
}
inline const IOdictionary& cfdemCloud::couplingProperties() const inline const IOdictionary& cfdemCloud::couplingProperties() const
{ {
return couplingProperties_; return couplingProperties_;
@ -311,6 +316,27 @@ inline label cfdemCloud::particleType(label index) const
} }
} }
inline bool cfdemCloud::getParticleAngVels() const
{
return getParticleAngVels_;
}
inline double ** cfdemCloud::particleAngVels() const
{
return particleAngVels_;
}
inline vector cfdemCloud::particleAngVel(label index) const
{
if(!getParticleAngVels_) return vector::zero;
else
{
vector angVel;
for(int i=0;i<3;i++) angVel[i] = particleAngVels_[index][i];
return angVel;
}
}
inline int cfdemCloud::numberOfParticles() const inline int cfdemCloud::numberOfParticles() const
{ {
return numberOfParticles_; return numberOfParticles_;

View File

@ -20,7 +20,9 @@ License
#include "cfdemCloudEnergy.H" #include "cfdemCloudEnergy.H"
#include "energyModel.H" #include "energyModel.H"
#include "massTransferModel.H"
#include "thermCondModel.H" #include "thermCondModel.H"
#include "diffCoeffModel.H"
#include "chemistryModel.H" #include "chemistryModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -40,9 +42,12 @@ cfdemCloudEnergy::cfdemCloudEnergy
: :
cfdemCloud(mesh), cfdemCloud(mesh),
energyModels_(couplingProperties_.lookup("energyModels")), energyModels_(couplingProperties_.lookup("energyModels")),
massTransferModels_(couplingProperties_.lookup("massTransferModels")),
implicitEnergyModel_(false), implicitEnergyModel_(false),
implicitMassTransferModel_(false),
chemistryModels_(couplingProperties_.lookup("chemistryModels")), chemistryModels_(couplingProperties_.lookup("chemistryModels")),
energyModel_(nrEnergyModels()), energyModel_(nrEnergyModels()),
massTransferModel_(nrMassTransferModels()),
thermCondModel_ thermCondModel_
( (
thermCondModel::New thermCondModel::New
@ -51,6 +56,14 @@ cfdemCloudEnergy::cfdemCloudEnergy
*this *this
) )
), ),
diffCoeffModel_
(
diffCoeffModel::New
(
couplingProperties_,
*this
)
),
chemistryModel_(nrChemistryModels()) chemistryModel_(nrChemistryModels())
{ {
forAll(energyModels_, modeli) forAll(energyModels_, modeli)
@ -67,6 +80,20 @@ cfdemCloudEnergy::cfdemCloudEnergy
); );
} }
forAll(massTransferModels_, modeli)
{
massTransferModel_.set
(
modeli,
massTransferModel::New
(
couplingProperties_,
*this,
massTransferModels_[modeli]
)
);
}
forAll(chemistryModels_, modeli) forAll(chemistryModels_, modeli)
{ {
chemistryModel_.set chemistryModel_.set
@ -100,6 +127,14 @@ void cfdemCloudEnergy::calcEnergyContributions()
} }
} }
void cfdemCloudEnergy::calcMassContributions()
{
forAll(massTransferModel_, modeli)
{
massTransferModel_[modeli].calcMassContribution();
}
}
void cfdemCloudEnergy::speciesExecute() void cfdemCloudEnergy::speciesExecute()
{ {
forAll(chemistryModel_, modeli) forAll(chemistryModel_, modeli)
@ -114,6 +149,11 @@ label cfdemCloudEnergy::nrEnergyModels() const
return energyModels_.size(); return energyModels_.size();
} }
label cfdemCloudEnergy::nrMassTransferModels() const
{
return massTransferModels_.size();
}
int cfdemCloudEnergy::nrChemistryModels() int cfdemCloudEnergy::nrChemistryModels()
{ {
return chemistryModels_.size(); return chemistryModels_.size();
@ -124,11 +164,22 @@ bool cfdemCloudEnergy::implicitEnergyModel() const
return implicitEnergyModel_; return implicitEnergyModel_;
} }
bool cfdemCloudEnergy::implicitMassTransferModel() const
{
return implicitMassTransferModel_;
}
const energyModel& cfdemCloudEnergy::energyM(int i) const energyModel& cfdemCloudEnergy::energyM(int i)
{ {
return energyModel_[i]; return energyModel_[i];
} }
const massTransferModel& cfdemCloudEnergy::massTransferM(int i)
{
return massTransferModel_[i];
}
const chemistryModel& cfdemCloudEnergy::chemistryM(int i) const chemistryModel& cfdemCloudEnergy::chemistryM(int i)
{ {
return chemistryModel_[i]; return chemistryModel_[i];
@ -139,6 +190,11 @@ const thermCondModel& cfdemCloudEnergy::thermCondM()
return thermCondModel_; return thermCondModel_;
} }
const diffCoeffModel& cfdemCloudEnergy::diffCoeffM()
{
return diffCoeffModel_;
}
void cfdemCloudEnergy::energyContributions(volScalarField& Qsource) void cfdemCloudEnergy::energyContributions(volScalarField& Qsource)
{ {
Qsource.primitiveFieldRef()=0.0; Qsource.primitiveFieldRef()=0.0;
@ -159,6 +215,26 @@ void cfdemCloudEnergy::energyCoefficients(volScalarField& Qcoeff)
} }
} }
void cfdemCloudEnergy::massContributions(volScalarField& Sm)
{
Sm.primitiveFieldRef()=0.0;
Sm.boundaryFieldRef()=0.0;
forAll(massTransferModel_, modeli)
{
massTransferM(modeli).addMassContribution(Sm);
}
}
void cfdemCloudEnergy::massCoefficients(volScalarField& Smi)
{
Smi.primitiveFieldRef()=0.0;
Smi.boundaryFieldRef()=0.0;
forAll(massTransferModel_, modeli)
{
massTransferM(modeli).addMassCoefficient(Smi);
}
}
bool cfdemCloudEnergy::evolve bool cfdemCloudEnergy::evolve
( (
volScalarField& alpha, volScalarField& alpha,
@ -184,6 +260,13 @@ bool cfdemCloudEnergy::evolve
if(verbose_) Info << "speciesExecute done" << endl; if(verbose_) Info << "speciesExecute done" << endl;
clockM().stop("speciesExecute"); clockM().stop("speciesExecute");
// calculate mass contributions
clockM().start(33,"calcMassContributions");
if(verbose_) Info << "- calcMassContributions" << endl;
calcMassContributions();
if(verbose_) Info << "calcMassContributions done." << endl;
clockM().stop("calcMassContributions");
return true; return true;
} }
return false; return false;
@ -196,6 +279,10 @@ void cfdemCloudEnergy::postFlow()
{ {
energyModel_[modeli].postFlow(); energyModel_[modeli].postFlow();
} }
forAll(massTransferModel_, modeli)
{
massTransferModel_[modeli].postFlow();
}
} }
void cfdemCloudEnergy::solve() void cfdemCloudEnergy::solve()
@ -204,6 +291,10 @@ void cfdemCloudEnergy::solve()
{ {
energyModel_[modeli].solve(); energyModel_[modeli].solve();
} }
forAll(massTransferModel_, modeli)
{
massTransferModel_[modeli].solve();
}
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -40,7 +40,9 @@ namespace Foam
{ {
// forward declarations // forward declarations
class energyModel; class energyModel;
class massTransferModel;
class thermCondModel; class thermCondModel;
class diffCoeffModel;
class chemistryModel; class chemistryModel;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class cfdemCloudEnergy Declaration Class cfdemCloudEnergy Declaration
@ -54,24 +56,36 @@ protected:
const wordList energyModels_; const wordList energyModels_;
const wordList massTransferModels_;
bool implicitEnergyModel_; bool implicitEnergyModel_;
bool implicitMassTransferModel_;
const wordList chemistryModels_; const wordList chemistryModels_;
PtrList<energyModel> energyModel_; PtrList<energyModel> energyModel_;
PtrList<massTransferModel> massTransferModel_;
autoPtr<thermCondModel> thermCondModel_; autoPtr<thermCondModel> thermCondModel_;
autoPtr<diffCoeffModel> diffCoeffModel_;
PtrList<chemistryModel> chemistryModel_; PtrList<chemistryModel> chemistryModel_;
void calcEnergyContributions(); void calcEnergyContributions();
void calcMassContributions();
void speciesExecute(); void speciesExecute();
public: public:
friend class energyModel; friend class energyModel;
friend class massTransferModel;
friend class chemistryModel; friend class chemistryModel;
// Constructors // Constructors
@ -92,24 +106,38 @@ public:
const energyModel& energyM(int); const energyModel& energyM(int);
const massTransferModel& massTransferM(int);
const thermCondModel& thermCondM(); const thermCondModel& thermCondM();
const diffCoeffModel& diffCoeffM();
const chemistryModel& chemistryM(int); const chemistryModel& chemistryM(int);
label nrEnergyModels() const; label nrEnergyModels() const;
label nrMassTransferModels() const;
bool implicitEnergyModel() const; bool implicitEnergyModel() const;
bool implicitMassTransferModel() const;
int nrChemistryModels(); int nrChemistryModels();
inline const wordList& energyModels() const; inline const wordList& energyModels() const;
inline const wordList& massTransferModels() const;
inline const wordList& chemistryModels() const; inline const wordList& chemistryModels() const;
void energyContributions(volScalarField&); void energyContributions(volScalarField&);
void energyCoefficients(volScalarField&); void energyCoefficients(volScalarField&);
void massContributions(volScalarField&);
void massCoefficients(volScalarField&);
bool evolve(volScalarField&,volVectorField&,volVectorField&); bool evolve(volScalarField&,volVectorField&,volVectorField&);
void postFlow(); void postFlow();

View File

@ -30,6 +30,11 @@ inline const wordList& cfdemCloudEnergy::energyModels() const
return energyModels_; return energyModels_;
} }
inline const wordList& cfdemCloudEnergy::massTransferModels() const
{
return massTransferModels_;
}
inline const wordList& cfdemCloudEnergy::chemistryModels() const inline const wordList& cfdemCloudEnergy::chemistryModels() const
{ {
return chemistryModels_; return chemistryModels_;

View File

@ -0,0 +1,108 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "SyamlalDiffCoeff.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(SyamlalDiffCoeff, 0);
addToRunTimeSelectionTable
(
diffCoeffModel,
SyamlalDiffCoeff,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
SyamlalDiffCoeff::SyamlalDiffCoeff
(
const dictionary& dict,
cfdemCloud& sm
)
:
diffCoeffModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_))
{
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
SyamlalDiffCoeff::~SyamlalDiffCoeff()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volScalarField> SyamlalDiffCoeff::diffCoeff() const
{
const volScalarField& D0Field_ = D0Field();
tmp<volScalarField> tvf
(
new volScalarField
(
IOobject
(
"tmpDiffCoeff",
voidfraction_.instance(),
voidfraction_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
voidfraction_.mesh(),
dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
)
);
volScalarField& svf = tvf.ref();
forAll(svf,cellI)
{
scalar D0 = D0Field_[cellI];
if (1-voidfraction_[cellI] < SMALL) svf[cellI] = D0;
else if (voidfraction_[cellI] < SMALL) svf[cellI] = 0.0;
else svf[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * D0;
}
return tvf;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,95 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
diffusion coefficient of species in fluid phase in presence of particles
according to Syamlal, M. and Gidaspow, M. AIChE Journal 31.1 (1985)
Class
SyamlalDiffCoeff
SourceFiles
SyamlalDiffCoeff.C
\*---------------------------------------------------------------------------*/
#ifndef SyamlalDiffCoeff_H
#define SyamlalDiffCoeff_H
#include "diffCoeffModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class SyamlalDiffCoeff Declaration
\*---------------------------------------------------------------------------*/
class SyamlalDiffCoeff
:
public diffCoeffModel
{
private:
dictionary propsDict_;
word voidfractionFieldName_;
const volScalarField& voidfraction_;
public:
//- Runtime type information
TypeName("SyamlalDiffCoeff");
// Constructors
//- Construct from components
SyamlalDiffCoeff
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~SyamlalDiffCoeff();
// Member Functions
tmp<volScalarField> diffCoeff() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "diffCoeffModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(diffCoeffModel, 0);
defineRunTimeSelectionTable(diffCoeffModel, dictionary);
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
diffCoeffModel::diffCoeffModel
(
const dictionary& dict,
cfdemCloud& sm
)
:
dict_(dict),
particleCloud_(sm),
transportProperties_
(
IOobject
(
"transportProperties",
sm.mesh().time().constant(),
sm.mesh(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
D0Field_
(
IOobject
(
"D0",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.)
)
{
// build constant fields for single phase case
if (!particleCloud_.multiphase())
{
D0Field_ = volScalarField
(
IOobject
(
"D0",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("D"))
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
diffCoeffModel::~diffCoeffModel()
{}
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
const volScalarField& diffCoeffModel::D0Field() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("D");
}
else
{
return D0Field_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,121 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
model for species diffusion coefficient in fluid phase in presence of particles
Class
diffCoeffModel
SourceFiles
diffCoeffModel.C
\*---------------------------------------------------------------------------*/
#ifndef diffCoeffModel_H
#define diffCoeffModel_H
#include "fvCFD.H"
#include "cfdemCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class diffCoeffModel Declaration
\*---------------------------------------------------------------------------*/
class diffCoeffModel
{
protected:
// Protected data
const dictionary& dict_;
cfdemCloud& particleCloud_;
IOdictionary transportProperties_;
mutable volScalarField D0Field_;
public:
//- Runtime type information
TypeName("diffCoeffModel");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
diffCoeffModel,
dictionary,
(
const dictionary& dict,
cfdemCloud& sm
),
(dict,sm)
);
// Constructors
//- Construct from components
diffCoeffModel
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
virtual ~diffCoeffModel();
// Selector
static autoPtr<diffCoeffModel> New
(
const dictionary& dict,
cfdemCloud& sm
);
// Member Functions
virtual tmp<volScalarField> diffCoeff() const = 0;
const volScalarField& D0Field() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,72 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "diffCoeffModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
autoPtr<diffCoeffModel> diffCoeffModel::New
(
const dictionary& dict,
cfdemCloud& sm
)
{
word diffCoeffModelType
(
dict.lookup("diffCoeffModel")
);
Info<< "Selecting diffCoeffModel "
<< diffCoeffModelType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(diffCoeffModelType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalError
<< "diffCoeffModel::New(const dictionary&, const spray&) : "
<< endl
<< " unknown diffCoeffModelType type "
<< diffCoeffModelType
<< ", constructor not in hash table" << endl << endl
<< " Valid diffCoeffModel types are :"
<< endl;
Info<< dictionaryConstructorTablePtr_->toc()
<< abort(FatalError);
}
return autoPtr<diffCoeffModel>(cstrIter()(dict,sm));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,89 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "noDiffCoeff.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(noDiffCoeff, 0);
addToRunTimeSelectionTable(diffCoeffModel, noDiffCoeff, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
noDiffCoeff::noDiffCoeff
(
const dictionary& dict,
cfdemCloud& sm
)
:
diffCoeffModel(dict,sm)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
noDiffCoeff::~noDiffCoeff()
{}
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
tmp<volScalarField> noDiffCoeff::diffCoeff() const
{
tmp<volScalarField> dcoeff
(
new volScalarField
(
IOobject
(
"fake1",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar
(
"zero",
dimensionSet(0,2,-1,0,0,0,0),
0.0
)
)
);
return dcoeff;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
Class
noDiffCoeff
SourceFiles
noDiffCoeff.C
\*---------------------------------------------------------------------------*/
#ifndef noDiffCoeff_H
#define noDiffCoeff_H
#include "diffCoeffModel.H"
#include "cfdemCloud.H"
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class noDrag Declaration
\*---------------------------------------------------------------------------*/
class noDiffCoeff
:
public diffCoeffModel
{
public:
//- Runtime type information
TypeName("off");
// Constructors
//- Construct from components
noDiffCoeff
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~noDiffCoeff();
// Member Functions
tmp<volScalarField> diffCoeff() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -55,15 +55,65 @@ energyModel::energyModel
IOobject::NO_WRITE IOobject::NO_WRITE
) )
), ),
kf0_ kf0Field_
( (
dimensionedScalar(transportProperties_.lookup("kf")).value() IOobject
(
"kf0",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
), ),
Cp_ sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.)
),
CpField_
( (
dimensionedScalar(transportProperties_.lookup("Cp")).value() IOobject
(
"Cp",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.)
) )
{} {
// build constant fields for single phase case
if (!particleCloud_.multiphase())
{
kf0Field_ = volScalarField
(
IOobject
(
"kf0",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("kf"))
);
CpField_ = volScalarField
(
IOobject
(
"Cp",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("Cp"))
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
@ -72,9 +122,28 @@ energyModel::~energyModel()
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
scalar energyModel::Cp() const const volScalarField& energyModel::kf0Field() const
{ {
return Cp_; if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("kf");
}
else
{
return kf0Field_;
}
}
const volScalarField& energyModel::CpField() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("Cp");
}
else
{
return CpField_;
}
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -46,9 +46,9 @@ protected:
IOdictionary transportProperties_; IOdictionary transportProperties_;
scalar kf0_; // fluid thermal conductivity [W/(m*K)] mutable volScalarField kf0Field_;
scalar Cp_; // specific heat capacity [W*s/(kg*K)] mutable volScalarField CpField_;
public: public:
@ -107,7 +107,9 @@ public:
virtual void solve() {} virtual void solve() {}
scalar Cp() const; const volScalarField& kf0Field() const;
const volScalarField& CpField() const;
virtual scalar aveTpart() const virtual scalar aveTpart() const
{ {

View File

@ -267,6 +267,9 @@ void heatTransferGunn::calcEnergyContribution()
const volScalarField mufField = particleCloud_.turbulence().nu()*rho_; const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
#endif #endif
const volScalarField& CpField_ = CpField();
const volScalarField& kf0Field_ = kf0Field();
if (typeCG_.size()>1 || typeCG_[0] > 1) if (typeCG_.size()>1 || typeCG_[0] > 1)
{ {
Info << "heatTransferGunn using scale = " << typeCG_ << endl; Info << "heatTransferGunn using scale = " << typeCG_ << endl;
@ -338,6 +341,9 @@ void heatTransferGunn::calcEnergyContribution()
ds = 2.*particleCloud_.radius(index); ds = 2.*particleCloud_.radius(index);
ds_scaled = ds/cg; ds_scaled = ds/cg;
scalar kf0 = kf0Field_[cellI];
scalar Cp = CpField_[cellI];
if (expNusselt_) if (expNusselt_)
{ {
Nup = NuField_[cellI]; Nup = NuField_[cellI];
@ -350,7 +356,7 @@ void heatTransferGunn::calcEnergyContribution()
magUr = mag(Ufluid - Us); magUr = mag(Ufluid - Us);
muf = mufField[cellI]; muf = mufField[cellI];
Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf; Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf;
Pr = max(SMALL, Cp_ * muf / kf0_); Pr = max(SMALL, Cp * muf / kf0);
Nup = Nusselt(voidfraction, Rep, Pr); Nup = Nusselt(voidfraction, Rep, Pr);
} }
Nup *= NusseltScalingFactor_; Nup *= NusseltScalingFactor_;
@ -358,7 +364,7 @@ void heatTransferGunn::calcEnergyContribution()
Tsum += partTemp_[index][0]; Tsum += partTemp_[index][0];
Nsum += 1.0; Nsum += 1.0;
scalar h = kf0_ * Nup / ds_scaled; scalar h = kf0 * Nup / ds_scaled;
scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere
// calc convective heat flux [W] // calc convective heat flux [W]
@ -380,8 +386,8 @@ void heatTransferGunn::calcEnergyContribution()
{ {
Pout << "partHeatFlux = " << partHeatFlux_[index][0] << endl; Pout << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
Pout << "magUr = " << magUr << endl; Pout << "magUr = " << magUr << endl;
Pout << "kf0 = " << kf0_ << endl; Pout << "kf0 = " << kf0 << endl;
Pout << "Cp = " << Cp_ << endl; Pout << "Cp = " << Cp << endl;
Pout << "rho = " << rho_[cellI] << endl; Pout << "rho = " << rho_[cellI] << endl;
Pout << "h = " << h << endl; Pout << "h = " << h << endl;
Pout << "ds = " << ds << endl; Pout << "ds = " << ds << endl;

View File

@ -266,6 +266,9 @@ void heatTransferRanzMarshall::calcEnergyContribution()
const volScalarField mufField = particleCloud_.turbulence().nu()*rho_; const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
#endif #endif
const volScalarField& CpField_ = CpField();
const volScalarField& kf0Field_ = kf0Field();
if (typeCG_.size()>1 || typeCG_[0] > 1) if (typeCG_.size()>1 || typeCG_[0] > 1)
{ {
Info << "heatTransferRanzMarshall using scale = " << typeCG_ << endl; Info << "heatTransferRanzMarshall using scale = " << typeCG_ << endl;
@ -338,6 +341,9 @@ void heatTransferRanzMarshall::calcEnergyContribution()
ds = 2.*particleCloud_.radius(index); ds = 2.*particleCloud_.radius(index);
ds_scaled = ds/cg; ds_scaled = ds/cg;
scalar kf0 = kf0Field_[cellI];
scalar Cp = CpField_[cellI];
if (expNusselt_) if (expNusselt_)
{ {
Nup = NuField_[cellI]; Nup = NuField_[cellI];
@ -350,7 +356,7 @@ void heatTransferRanzMarshall::calcEnergyContribution()
magUr = mag(Ufluid - Us); magUr = mag(Ufluid - Us);
muf = mufField[cellI]; muf = mufField[cellI];
Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf; Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf;
Pr = max(SMALL, Cp_ * muf / kf0_); Pr = max(SMALL, Cp * muf / kf0);
Nup = Nusselt(voidfraction, Rep, Pr); Nup = Nusselt(voidfraction, Rep, Pr);
} }
Nup *= NusseltScalingFactor_; Nup *= NusseltScalingFactor_;
@ -358,7 +364,7 @@ void heatTransferRanzMarshall::calcEnergyContribution()
Tsum += partTemp_[index][0]; Tsum += partTemp_[index][0];
Nsum += 1.0; Nsum += 1.0;
scalar h = kf0_ * Nup / ds_scaled; scalar h = kf0 * Nup / ds_scaled;
scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere
// calc convective heat flux [W] // calc convective heat flux [W]

View File

@ -0,0 +1,376 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "wallHeatTransferYagi.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(wallHeatTransferYagi, 0);
addToRunTimeSelectionTable(energyModel, wallHeatTransferYagi, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
wallHeatTransferYagi::wallHeatTransferYagi
(
const dictionary& dict,
cfdemCloudEnergy& sm
)
:
energyModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
QWallFluidName_(propsDict_.lookupOrDefault<word>("QWallFluidName","QWallFluid")),
QWallFluid_
( IOobject
(
QWallFluidName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
),
QWallFluidCoeffName_(propsDict_.lookupOrDefault<word>("QWallFluidCoeffName","QWallFluidCoeff")),
QWallFluidCoeff_
( IOobject
(
QWallFluidCoeffName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
),
wallTempName_(propsDict_.lookup("wallTempName")),
wallTemp_
( IOobject
(
wallTempName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0)
),
dpField_
( IOobject
(
"dpField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,1,0,0,0,0,0), 0.0)
),
GField_
( IOobject
(
"GField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(1,-2,-1,0,0,0,0), vector(0.0, 0.0, 0.0))
),
magGField_
( IOobject
(
"magGField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-2,-1,0,0,0,0), 0.0)
),
ReField_
( IOobject
(
"ReField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
),
PrField_
( IOobject
(
"PrField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
),
tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
voidfractionMax_(readScalar(propsDict_.lookup("voidfractionMax"))),
maxSource_(1e30),
velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
dpArrayRegName_(typeName + "dpArray")
{
particleCloud_.registerParticleProperty<double**>(dpArrayRegName_,1);
if (propsDict_.found("maxSource"))
{
maxSource_=readScalar(propsDict_.lookup ("maxSource"));
Info << "limiting wall source field to: " << maxSource_ << endl;
}
if (!implicit_)
{
QWallFluidCoeff_.writeOpt() = IOobject::NO_WRITE;
}
if (verbose_)
{
dpField_.writeOpt() = IOobject::AUTO_WRITE;
GField_.writeOpt() = IOobject::AUTO_WRITE;
magGField_.writeOpt() = IOobject::AUTO_WRITE;
ReField_.writeOpt() = IOobject::AUTO_WRITE;
PrField_.writeOpt() = IOobject::AUTO_WRITE;
dpField_.write();
GField_.write();
magGField_.write();
ReField_.write();
PrField_.write();
}
// currently it is detected if field was auto generated or defined
// improvement would be changing the type here automatically
forAll(wallTemp_.boundaryField(),patchI)
if(wallTemp_.boundaryField()[patchI].type()=="calculated")
FatalError <<"Scalar field:"<< wallTemp_.name() << " must be defined.\n" << abort(FatalError);
wallTemp_.writeOpt() = IOobject::AUTO_WRITE;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
wallHeatTransferYagi::~wallHeatTransferYagi()
{
}
// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
void wallHeatTransferYagi::calcEnergyContribution()
{
double**& dpArray_ = particleCloud_.getParticlePropertyRef<double**>(dpArrayRegName_);
// reset Scalar field
QWallFluid_.primitiveFieldRef() = 0.0;
if (implicit_)
QWallFluidCoeff_.primitiveFieldRef() = 0.0;
#ifdef compre
const volScalarField mufField = particleCloud_.turbulence().mu();
#else
const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
#endif
const volScalarField& CpField_ = CpField();
const volScalarField& kf0Field_ = kf0Field();
// calculate mean dp
for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
{
label cellI = particleCloud_.cellIDs()[index][0];
if(cellI >= 0)
{
scalar ds = 2.*particleCloud_.radius(index);
dpArray_[index][0] = ds;
}
}
// calculate mean dp field
particleCloud_.averagingM().resetWeightFields();
particleCloud_.averagingM().setScalarAverage
(
dpField_,
dpArray_,
particleCloud_.particleWeights(),
particleCloud_.averagingM().UsWeightField(),
NULL
);
// calculate G field (superficial mass velocity)
GField_ = U_*voidfraction_*rho_;
magGField_ = mag(GField_);
// calculate Re field
ReField_ = dpField_ * magGField_ / mufField;
// calculate Pr field
PrField_ = CpField_ * mufField / kf0Field_;
const fvPatchList& patches = U_.mesh().boundary();
// calculate flux
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (wallTemp_.boundaryField().types()[patchi] == "fixedValue")
{
if(tempField_.boundaryField().types()[patchi] == "zeroGradient")
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
// calculate H
scalar JH;
if (voidfraction_[faceCelli]<=voidfractionMax_)
JH = 0.20 * pow(ReField_[faceCelli]+SMALL,-0.20); // Yagi eq. 6a
else
JH = 0;
scalar h = JH * CpField_[faceCelli] * magGField_[faceCelli] / (pow(PrField_[faceCelli],0.666) + SMALL);
// get delta T (wall-fluid)
scalar Twall = wallTemp_.boundaryField()[patchi][facei];
scalar Tfluid = tempField_[faceCelli];
// get area
scalar area = curPatch.magSf()[facei];
// calculate heat flux
heatFlux(faceCelli, h, area, Twall, Tfluid);
if(verbose_ && facei >=0 && facei <2)
{
scalar deltaT = Twall - Tfluid;
Info << "####################" << endl;
Info << "cellID: " << faceCelli << endl;
Info << "kf: " << kf0Field_[faceCelli] << " J/msK" << endl;
Info << "Cp: " << CpField_[faceCelli] << " J/kgK" << endl;
Info << "ro: " << rho_[faceCelli] << " kg/m3" << endl;
Info << "mu: " << mufField[faceCelli] << " Pa s" << endl;
Info << "dp: " << dpField_[faceCelli] << " m" << endl;
Info << "ep: " << voidfraction_[faceCelli] << endl;
Info << "U : " << U_[faceCelli] << " m/s" << endl;
Info << "G : " << GField_[faceCelli] << " kg/m2s" << endl;
Info << "mG: " << magGField_[faceCelli] << " kg/m2s" << endl;
Info << "Re: " << ReField_[faceCelli] << endl;
Info << "Pr: " << PrField_[faceCelli] << endl;
Info << "JH: " << JH << endl;
Info << "h : " << h << " J/m2sK" << endl;
Info << "Tw: " << Twall << " K" << endl;
Info << "Tf: " << Tfluid << " K" << endl;
Info << "dT: " << deltaT << " K" << endl;
Info << "q : " << h*deltaT << " J/m2s" << endl;
Info << "A : " << area << " m2" << endl;
Info << "Q : " << h*deltaT*area << " J/s" << endl;
}
}
}
else
{
FatalError << "wallHeatTransferYagi requires zeroGradient BC for temperature field" << endl;
}
}
}
QWallFluid_.primitiveFieldRef() /= QWallFluid_.mesh().V();
if(implicit_)
QWallFluidCoeff_.primitiveFieldRef() /= QWallFluidCoeff_.mesh().V();
// limit source term in explicit treatment
if(!implicit_)
{
forAll(QWallFluid_,cellI)
{
scalar EuFieldInCell = QWallFluid_[cellI];
if(mag(EuFieldInCell) > maxSource_ )
{
Pout << "limiting source term" << endl ;
QWallFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
}
}
}
QWallFluid_.correctBoundaryConditions();
}
void wallHeatTransferYagi::addEnergyContribution(volScalarField& Qsource) const
{
Qsource += QWallFluid_;
}
void wallHeatTransferYagi::addEnergyCoefficient(volScalarField& Qcoeff) const
{
if(implicit_)
{
Qcoeff += QWallFluidCoeff_;
}
}
void wallHeatTransferYagi::heatFlux(label faceCelli, scalar h, scalar area, scalar Twall, scalar Tfluid)
{
if(!implicit_)
{
QWallFluid_[faceCelli] += h * area * (Twall - Tfluid);
}
else
{
QWallFluid_[faceCelli] += h * area * Twall;
QWallFluidCoeff_[faceCelli] -= h * area;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
Correlation for wall-to-bed heat transfer according to
Yagi, S. A.I.Ch.E. Journal 5.1 (1959)
\*---------------------------------------------------------------------------*/
#ifndef wallHeatTransferYagi_H
#define wallHeatTransferYagi_H
#include "fvCFD.H"
#include "cfdemCloudEnergy.H"
#include "energyModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class wallHeatTransferYagi Declaration
\*---------------------------------------------------------------------------*/
class wallHeatTransferYagi
:
public energyModel
{
protected:
dictionary propsDict_;
bool verbose_;
bool implicit_;
word QWallFluidName_;
volScalarField QWallFluid_;
word QWallFluidCoeffName_;
volScalarField QWallFluidCoeff_;
word wallTempName_;
volScalarField wallTemp_;
volScalarField dpField_;
volVectorField GField_;
volScalarField magGField_;
volScalarField ReField_;
volScalarField PrField_;
word tempFieldName_;
const volScalarField& tempField_;
word voidfractionFieldName_;
const volScalarField& voidfraction_;
scalar voidfractionMax_;
scalar maxSource_; // max (limited) value of src field
word velFieldName_;
const volVectorField& U_;
word densityFieldName_;
const volScalarField& rho_;
const word dpArrayRegName_;
virtual void heatFlux(label, scalar, scalar, scalar, scalar);
public:
//- Runtime type information
TypeName("wallHeatTransferYagi");
// Constructors
//- Construct from components
wallHeatTransferYagi
(
const dictionary& dict,
cfdemCloudEnergy& sm
);
// Destructor
virtual ~wallHeatTransferYagi();
// Member Functions
void addEnergyContribution(volScalarField&) const;
void addEnergyCoefficient(volScalarField&) const;
void calcEnergyContribution();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -209,6 +209,7 @@ void BeetstraDrag::setForce() const
Ufluid = UInterpolator_.interpolate(position,cellI); Ufluid = UInterpolator_.interpolate(position,cellI);
//Ensure interpolated void fraction to be meaningful //Ensure interpolated void fraction to be meaningful
// Info << " --> voidfraction: " << voidfraction << endl; // Info << " --> voidfraction: " << voidfraction << endl;
if (voidfraction > 1.00) voidfraction = 1.0; if (voidfraction > 1.00) voidfraction = 1.0;
if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_; if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
} }
@ -280,6 +281,7 @@ void BeetstraDrag::setForce() const
{ {
Pout << "cellI = " << cellI << endl; Pout << "cellI = " << cellI << endl;
Pout << "index = " << index << endl; Pout << "index = " << index << endl;
Pout << "Ufluid = " << Ufluid << endl;
Pout << "Us = " << Us << endl; Pout << "Us = " << Us << endl;
Pout << "Ur = " << Ur << endl; Pout << "Ur = " << Ur << endl;
Pout << "ds = " << ds << endl; Pout << "ds = " << ds << endl;

View File

@ -102,8 +102,6 @@ protected:
double h(double) const; double h(double) const;
public: public:
//- Runtime type information //- Runtime type information

View File

@ -81,9 +81,29 @@ MeiLift::MeiLift
propsDict_(dict.subDict(typeName + "Props")), propsDict_(dict.subDict(typeName + "Props")),
velFieldName_(propsDict_.lookup("velFieldName")), velFieldName_(propsDict_.lookup("velFieldName")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)), U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
useSecondOrderTerms_(false) useShearInduced_(propsDict_.lookupOrDefault<bool>("useShearInduced",true)),
useSpinInduced_(propsDict_.lookupOrDefault<bool>("useSpinInduced",false)),
combineShearSpin_(propsDict_.lookupOrDefault<bool>("combineShearSpin",false))
{ {
if (propsDict_.found("useSecondOrderTerms")) useSecondOrderTerms_=true; // read switches
if(useShearInduced_)
Info << "Lift model: including shear-induced term.\n";
if(useSpinInduced_)
{
Info << "Lift model: including spin-induced term.\n";
Info << "Make sure to use a rolling friction model in LIGGGHTS!\n";
if(!dict.lookupOrDefault<bool>("getParticleAngVels",false))
FatalError << "Lift model: useSpinInduced=true requires getParticleAngVels=true in couplingProperties" << abort(FatalError);
}
if(combineShearSpin_)
{
Info << "Lift model: combining shear- and spin-terms by assuming equilibrium spin-rate.\n";
if(!useShearInduced_ || !useSpinInduced_)
FatalError << "Shear- and spin-induced lift must be activated in order to combine." << abort(FatalError);
}
// init force sub model // init force sub model
setForceSubModels(propsDict_); setForceSubModels(propsDict_);
@ -100,10 +120,15 @@ MeiLift::MeiLift
particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().initialize(typeName, typeName+".logDat");
particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force
particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug
particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug particleCloud_.probeM().vectorFields_.append("vorticity");
particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug particleCloud_.probeM().vectorFields_.append("Ang_velocity");
particleCloud_.probeM().scalarFields_.append("Rew"); //other are debug particleCloud_.probeM().scalarFields_.append("Rep");
particleCloud_.probeM().scalarFields_.append("J_star"); //other are debug particleCloud_.probeM().scalarFields_.append("Rew");
particleCloud_.probeM().scalarFields_.append("J*");
particleCloud_.probeM().scalarFields_.append("Cl(shear)");
particleCloud_.probeM().scalarFields_.append("Cl*(spin)");
particleCloud_.probeM().scalarFields_.append("Omega_eq");
particleCloud_.probeM().scalarFields_.append("Cl(combined)");
particleCloud_.probeM().writeHeader(); particleCloud_.probeM().writeHeader();
} }
@ -122,27 +147,44 @@ void MeiLift::setForce() const
const volScalarField& nufField = forceSubM(0).nuField(); const volScalarField& nufField = forceSubM(0).nuField();
const volScalarField& rhoField = forceSubM(0).rhoField(); const volScalarField& rhoField = forceSubM(0).rhoField();
// vectors
vector position(0,0,0); vector position(0,0,0);
vector lift(0,0,0); vector lift(0,0,0);
vector Us(0,0,0); vector Us(0,0,0);
vector Ur(0,0,0); vector Ur(0,0,0);
vector Omega(0,0,0);
vector vorticity(0,0,0);
// properties
scalar magUr(0); scalar magUr(0);
scalar magVorticity(0); scalar magVorticity(0);
scalar magOmega(0);
scalar ds(0); scalar ds(0);
scalar nuf(0); scalar nuf(0);
scalar rho(0); scalar rho(0);
scalar voidfraction(1); scalar voidfraction(1);
// dimensionless groups
scalar Rep(0); scalar Rep(0);
scalar Rew(0); scalar Rew(0);
scalar Cl(0);
scalar Cl_star(0); // shear induced
scalar J_star(0);
scalar Omega_eq(0);
scalar alphaStar(0);
scalar epsilonSqr(0.0);
scalar epsilon(0);
scalar omega_star(0); scalar omega_star(0);
vector vorticity(0,0,0); scalar Clshear(0);
scalar J_star(0);
scalar alphaStar(0);
scalar epsilonSqr(0);
scalar epsilon(0);
// spin induced
scalar Omega_star(0);
scalar Clspin_star(0);
// shear-spin combination
scalar Omega_eq(0);
scalar Clcombined(0);
volVectorField vorticityField = fvc::curl(U_); volVectorField vorticityField = fvc::curl(U_);
interpolationCellPoint<vector> UInterpolator_(U_); interpolationCellPoint<vector> UInterpolator_(U_);
@ -152,14 +194,14 @@ void MeiLift::setForce() const
for (int index = 0; index < particleCloud_.numberOfParticles(); ++index) for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
{ {
//if(mask[index][0])
//{
lift = vector::zero; lift = vector::zero;
label cellI = particleCloud_.cellIDs()[index][0]; label cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1) // particle Found if (cellI > -1) // particle Found
{ {
// properties
Us = particleCloud_.velocity(index); Us = particleCloud_.velocity(index);
Omega = particleCloud_.particleAngVel(index);
if (forceSubM(0).interpolation()) if (forceSubM(0).interpolation())
{ {
@ -173,24 +215,24 @@ void MeiLift::setForce() const
vorticity = vorticityField[cellI]; vorticity = vorticityField[cellI];
} }
magUr = mag(Ur);
magVorticity = mag(vorticity);
if (magUr > 0 && magVorticity > 0)
{
ds = 2. * particleCloud_.radius(index); ds = 2. * particleCloud_.radius(index);
nuf = nufField[cellI]; nuf = nufField[cellI];
rho = rhoField[cellI]; rho = rhoField[cellI];
// calc dimensionless properties magUr = mag(Ur);
Rep = ds*magUr/nuf; Rep = ds*magUr/nuf;
Rew = magVorticity*ds*ds/nuf;
// shear-induced lift
if (useShearInduced_)
{
magVorticity = mag(vorticity);
Rew = magVorticity*ds*ds/nuf;
omega_star = magVorticity * ds / magUr; omega_star = magVorticity * ds / magUr;
alphaStar = 0.5 * omega_star; alphaStar = 0.5 * omega_star;
epsilonSqr = omega_star / Rep; epsilonSqr = omega_star / Rep;
epsilon = sqrt(epsilonSqr); epsilon = sqrt(epsilonSqr);
//Basic model for the correction to the Saffman lift //Basic model for the correction to the Saffman lift
//McLaughlin (1991), Mei (1992), Loth and Dorgan (2009) //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009)
//J_star = 0.443 * J //J_star = 0.443 * J
@ -212,35 +254,76 @@ void MeiLift::setForce() const
* (1.0 + tanh(2.5 * (log10(epsilon) + 0.191))) * (1.0 + tanh(2.5 * (log10(epsilon) + 0.191)))
* (0.667 + tanh(6.0 * ( epsilon - 0.32 ))); * (0.667 + tanh(6.0 * ( epsilon - 0.32 )));
} }
//Loth and Dorgan (2009), Eq (31): Saffman lift-coefficient: ClSaff = 12.92 / pi * epsilon ~ 4.11 * epsilon
//Loth and Dorgan (2009), Eq (32)
Cl = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model
//Second order terms given by Loth and Dorgan (2009) //Loth and Dorgan (2009), Eq (31), Eq (32)
if (useSecondOrderTerms_) Clshear = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model
{
scalar sqrtRep = sqrt(Rep);
//Loth and Dorgan (2009), Eq (34)
Cl_star = 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (alphaStar - 2.0)))) * tanh(0.18 * sqrtRep);
//Loth and Dorgan (2009), Eq (38)
Omega_eq = alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrtRep - 0.001 * Rep);
//Loth and Dorgan (2009), Eq (39)
Cl += Omega_eq * Cl_star;
} }
if (useSpinInduced_)
{
magOmega = mag(Omega);
Omega_star = magOmega * ds / magUr;
//Loth and Dorgan (2009), Eq (34)
Clspin_star = 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (Omega_star - 2.0)))) * tanh(0.18 * sqrt(Rep));
}
if (combineShearSpin_)
{
//Loth and Dorgan (2009), Eq (38)
Omega_eq = alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrt(Rep) - 0.001 * Rep);
//Loth and Dorgan (2009), Eq (39)
Clcombined = Clshear + Clspin_star * Omega_eq;
if (magUr>0.0 && magVorticity>0.0)
{
//Loth and Dorgan (2009), Eq (27) //Loth and Dorgan (2009), Eq (27)
// please note: in Loth and Dorgan Vrel = Vp-Uf, here Urel = Uf-Vp. Hence the reversed cross products.
lift = 0.125 * constant::mathematical::pi lift = 0.125 * constant::mathematical::pi
* rho * rho
* Cl * Clcombined
* magUr * Ur ^ vorticity / magVorticity * magUr * magUr
* (Ur ^ vorticity) / mag(Ur ^ vorticity) // force direction
* ds * ds; * ds * ds;
}
}
else
{
//Loth and Dorgan (2009), Eq (36)
// please note: in Loth and Dorgan Vrel = Vp-Uf, here Urel = Uf-Vp. Hence the reversed cross products.
if (useShearInduced_)
{
if (magUr>0.0 && magVorticity>0.0)
{
//Loth and Dorgan (2009), Eq (27)
lift += 0.125 * constant::mathematical::pi
* rho
* Clshear
* magUr * magUr
* (Ur ^ vorticity) / mag(Ur ^ vorticity) // force direction
* ds * ds;
}
}
if (useSpinInduced_)
{
if (magUr>0.0 && magOmega>0.0)
{
//Loth and Dorgan (2009), Eq (33)
lift += 0.125 * constant::mathematical::pi
* rho
* Clspin_star
* (Ur ^ Omega)
* ds * ds * ds;
}
}
}
if (modelType_ == "B") if (modelType_ == "B")
{ {
voidfraction = particleCloud_.voidfraction(index); voidfraction = particleCloud_.voidfraction(index);
lift /= voidfraction; lift /= voidfraction;
} }
}
//********************************** //**********************************
//SAMPLING AND VERBOSE OUTOUT //SAMPLING AND VERBOSE OUTOUT
@ -250,6 +333,7 @@ void MeiLift::setForce() const
Pout << "Us = " << Us << endl; Pout << "Us = " << Us << endl;
Pout << "Ur = " << Ur << endl; Pout << "Ur = " << Ur << endl;
Pout << "vorticity = " << vorticity << endl; Pout << "vorticity = " << vorticity << endl;
Pout << "Omega = " << Omega << endl;
Pout << "ds = " << ds << endl; Pout << "ds = " << ds << endl;
Pout << "rho = " << rho << endl; Pout << "rho = " << rho << endl;
Pout << "nuf = " << nuf << endl; Pout << "nuf = " << nuf << endl;
@ -258,6 +342,10 @@ void MeiLift::setForce() const
Pout << "alphaStar = " << alphaStar << endl; Pout << "alphaStar = " << alphaStar << endl;
Pout << "epsilon = " << epsilon << endl; Pout << "epsilon = " << epsilon << endl;
Pout << "J_star = " << J_star << endl; Pout << "J_star = " << J_star << endl;
Pout << "Omega_eq = " << Omega_eq << endl;
Pout << "Clshear = " << Clshear<< endl;
Pout << "Clspin_star = " << Clspin_star << endl;
Pout << "Clcombined = " << Clcombined << endl;
Pout << "lift = " << lift << endl; Pout << "lift = " << lift << endl;
} }
@ -268,9 +356,14 @@ void MeiLift::setForce() const
vValues.append(lift); //first entry must the be the force vValues.append(lift); //first entry must the be the force
vValues.append(Ur); vValues.append(Ur);
vValues.append(vorticity); vValues.append(vorticity);
vValues.append(Omega);
sValues.append(Rep); sValues.append(Rep);
sValues.append(Rew); sValues.append(Rew);
sValues.append(J_star); sValues.append(J_star);
sValues.append(Clshear);
sValues.append(Clspin_star);
sValues.append(Omega_eq);
sValues.append(Clcombined);
particleCloud_.probeM().writeProbe(index, sValues, vValues); particleCloud_.probeM().writeProbe(index, sValues, vValues);
} }
// END OF SAMPLING AND VERBOSE OUTOUT // END OF SAMPLING AND VERBOSE OUTOUT
@ -279,7 +372,6 @@ void MeiLift::setForce() const
} }
// write particle based data to global array // write particle based data to global array
forceSubM(0).partToArray(index,lift,vector::zero); forceSubM(0).partToArray(index,lift,vector::zero);
//}
} }
} }

View File

@ -83,7 +83,11 @@ private:
const volVectorField& U_; const volVectorField& U_;
bool useSecondOrderTerms_; bool useShearInduced_;
bool useSpinInduced_;
bool combineShearSpin_;
public: public:

View File

@ -0,0 +1,624 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "ParmarBassetForce.H"
#include "addToRunTimeSelectionTable.H"
#include "smoothingModel.H"
#include "constDiffSmoothing.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
#define NOTONCPU 9999
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(ParmarBassetForce, 0);
addToRunTimeSelectionTable
(
forceModel,
ParmarBassetForce,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
ParmarBassetForce::ParmarBassetForce
(
const dictionary& dict,
cfdemCloud& sm
)
:
forceModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
velFieldName_(propsDict_.lookup("velFieldName")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
UsFieldName_(propsDict_.lookup("granVelFieldName")),
Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
nInt_(readLabel(propsDict_.lookup("nIntegral"))),
discOrder_(readLabel(propsDict_.lookup("discretisationOrder"))),
nHist_(nInt_+discOrder_+1),
FHistSize_(2*discOrder_+1),
ddtUrelHist_(nHist_,NULL), // UrelHist_[ndt in past][particle ID][dim]
rHist_(nHist_,NULL), // rHist_[ndt in past][particle ID][0]
FHist_(2,List<double**>(FHistSize_,NULL)), // FHist_[k={1,2}-1][ndt in past][particle ID][dim]
gH0_(NULL),
tRef_(NULL),
mRef_(NULL),
lRef_(NULL),
Urel_
( IOobject
(
"Urel",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector::zero)
),
ddtUrel_
( IOobject
(
"ddtUrel",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(0,1,-2,0,0,0,0), vector::zero)
),
smoothingModel_
(
smoothingModel::New
(
propsDict_,
sm
)
)
{
// init force sub model
setForceSubModels(propsDict_);
// define switches which can be read from dict
forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
forceSubM(0).readSwitches();
//Extra switches/settings
particleCloud_.checkCG(true);
if (discOrder_ < 1 || discOrder_ > 2)
FatalError << "Parmar Basset Force: Discretisation order > 2 not implemented!" << abort(FatalError);
//Append the field names to be probed
particleCloud_.probeM().initialize(typeName, typeName+".logDat");
particleCloud_.probeM().vectorFields_.append("ParmarBassetForce"); //first entry must the be the force
particleCloud_.probeM().vectorFields_.append("Urel");
particleCloud_.probeM().vectorFields_.append("ddtUrel");
//
particleCloud_.probeM().vectorFields_.append("UrelNoSmooth");
particleCloud_.probeM().vectorFields_.append("ddtUrelNoSmooth");
//
particleCloud_.probeM().vectorFields_.append("Fshort");
particleCloud_.probeM().vectorFields_.append("Flong1");
particleCloud_.probeM().vectorFields_.append("Flong2");
particleCloud_.probeM().scalarFields_.append("ReRef");
particleCloud_.probeM().scalarFields_.append("tRef");
particleCloud_.probeM().scalarFields_.append("mRef");
particleCloud_.probeM().scalarFields_.append("lRef");
particleCloud_.probeM().scalarFields_.append("r");
particleCloud_.probeM().scalarFields_.append("t0");
particleCloud_.probeM().scalarFields_.append("nKnown");
particleCloud_.probeM().writeHeader();
Urel_ = Us_ - U_;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
ParmarBassetForce::~ParmarBassetForce()
{
particleCloud_.dataExchangeM().destroy(gH0_, 1);
particleCloud_.dataExchangeM().destroy(tRef_, 1);
particleCloud_.dataExchangeM().destroy(mRef_, 1);
particleCloud_.dataExchangeM().destroy(lRef_, 1);
for (int i=0; i<nHist_; i++)
{
particleCloud_.dataExchangeM().destroy(ddtUrelHist_[i],3);
particleCloud_.dataExchangeM().destroy(rHist_ [i],1);
}
for (int i=0; i<FHistSize_; i++)
for (int k=0; k<2; k++)
particleCloud_.dataExchangeM().destroy(FHist_[k][i],3);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void ParmarBassetForce::setForce() const
{
// allocate arrays
if(particleCloud_.numberOfParticlesChanged())
reAllocArrays();
vector position(0,0,0);
vector Urel(0,0,0);
vector ddtUrel(0,0,0);
scalar t0min = 0.00;
const volScalarField& nufField = forceSubM(0).nuField();
const volScalarField& rhoField = forceSubM(0).rhoField();
double c[2][4][3] = {{{0.004045,0.011788,-0.851409},
{0.184350,0.009672,-0.425408},
{0.000746,0.014390,-1.262288},
{0.021708,0.016611,-0.867352}},
{{0.536783,0.005453,-1.430163},
{1.683212,0.003598,-0.773872},
{0.316350,0.004367,-1.423649},
{0.420423,0.000677,-0.670551}}
};
double chi[2][4][2] = {{{1.464787, 0.422501},
{0.785093, -0.078301},
{0.585753, -0.118737},
{-0.077407,-0.016665}},
{{1.001435, -0.003566},
{0.415345, -0.136697},
{0.536055, -0.205041},
{-0.032671,-0.021047}}
};
#include "setupProbeModel.H"
Urel_ = Us_ - U_;
ddtUrel_ = fvc::ddt(Us_) - fvc::ddt(U_) - (Us_ & fvc::grad(U_));
//
volVectorField UrelNoSmooth_ = Urel_;
volVectorField ddtUrelNoSmooth_ = ddtUrel_;
//
smoothingM().smoothen(Urel_);
smoothingM().smoothen(ddtUrel_);
interpolationCellPoint<vector> UrelInterpolator_(Urel_);
interpolationCellPoint<vector> ddtUrelInterpolator_(ddtUrel_);
//
interpolationCellPoint<vector> UrelNoSmoothInterpolator_(UrelNoSmooth_);
interpolationCellPoint<vector> ddtUrelNoSmoothInterpolator_(ddtUrelNoSmooth_);
//
for(int index = 0;index < particleCloud_.numberOfParticles(); index++)
{
vector ParmarBassetForce(0,0,0);
vector Fshort(0,0,0);
label cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1) // particle Found
{
// particle position (m)
position = particleCloud_.position(index);
// relative velocity (m/s)
if(forceSubM(0).interpolation())
Urel = UrelInterpolator_.interpolate(position,cellI);
else
Urel = Urel_[cellI];
// acceleration (m/s2)
if(forceSubM(0).interpolation())
ddtUrel = ddtUrelInterpolator_.interpolate(position,cellI);
else
ddtUrel = ddtUrel_[cellI];
//********* set/get reference point *********//
scalar rs = particleCloud_.radius(index);
scalar Re = 2.*rs*mag(Urel)/nufField[cellI];
scalar gH = (0.75 + 0.105*Re)/(Re+SMALL); // Eq. 3.3
scalar r;
if (gH0_[index][0]!=NOTONCPU)
{
r = pow(gH0_[index][0]/gH,1.5); // Eq. 3.4
if (r<0.25 || r>2.0)
{
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
}
if (gH0_[index][0]==NOTONCPU)
{
// calculate new reference values
scalar tRef = (rs*rs) / nufField[cellI] * (gH*gH) * 4.3354; // (256/pi)^(1/3) = 4.3354 (Eq 3.1)
scalar Vs = rs*rs*rs*M_PI*4/3;
scalar mRef = Vs*rhoField[cellI] * gH * 5.2863; // 9/(2*sqrt(pi))*(256/pi)^(1/6) = 5.2863 (Eq. 3.2)
gH0_[index][0] = gH;
tRef_[index][0] = tRef;
mRef_[index][0] = mRef;
lRef_[index][0] = rs;
r = 1;
}
scalar mps = lRef_[index][0]/tRef_[index][0]; // m/s
scalar mpss = mps/tRef_[index][0]; // m/s^2
scalar newton = mpss*mRef_[index][0]; // kg*m/s^2
scalar dt = U_.mesh().time().deltaT().value() / tRef_[index][0]; // dim.less
scalar t0 = nInt_*dt; // dim.less
//********* update histories *********//
// non-dimensionlise
Urel /= mps;
ddtUrel /= mpss;
// update ddtUrel and r history
update_ddtUrelHist(ddtUrel,index); // add current dim.less ddtUrel to history
update_rHist(r,index); // add current r to history
// warning and reset for too small t0
if (t0<t0min)
{
Pout << "ParmarBassetForce WARNING: t0 = " << t0 << " at ID = " << index << endl;
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
// check length of known history
int nKnown = 0;
for (int j=0; j<nHist_; j++) // loop over past times
{
if (ddtUrelHist_[j][index][0] == NOTONCPU)
break;
else
nKnown++;
}
//********* short term force computing (full integral) *********//
// number of values to be used for short term force computing
int nShort = min(nKnown,nInt_+1);
// int_0^dt K(r,xi) dxi * ddtU(t) dxi (singularity treated by assuming constant acceleration)
if (nShort>0)
{
for (int i=0; i<3; i++) // loop over dimensions
Fshort[i] = -calculateK0(r,dt) * ddtUrelHist_[0][index][i];
}
// int_dt^t0 K(r,xi) * ddtU(t-xi) dxi (trapezoid rule)
if (nShort>2)
{
for (int j=1; j<nShort; j++)
{
scalar xi = j*dt;
scalar K = pow((pow(xi,.25) + rHist_[j][index][0]*xi),-2.); // Eq. 3.4
for (int i=0; i<3; i++) // loop over dimensions
Fshort[i] -= trapWeight(j,nShort) * K * ddtUrelHist_[j][index][i] * dt;
}
}
//********* long term force computing (differential form) *********//
// update F1, F2 history
update_FHist(vector::zero,vector::zero,index);
// initialise ddtUrel(t0) and Flong(:) as 0 and r(t0) as 1
if (nKnown == nInt_)
{
for (int j=nInt_; j<nHist_; j++) // loop over past times
{
rHist_[j][index][0] = 1.;
for (int i=0; i<3; i++) // loop over dimensions
ddtUrelHist_[j][index][i] = 0.0;
}
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=0; j<FHistSize_; j++) // loop over past times
for (int i=0; i<3; i++) // loop over dimensions
FHist_[k][j][index][i] = 0.0;
nKnown = nHist_;
}
// solve ODEs
if (nKnown == nHist_)
{
for (int k=0; k<2; k++) // loop over F1, F2
{
//calculate coefficients
double C[4];
calculateCoeffs(k,t0,rHist_[nInt_][index][0],c,chi,C);
// solve Eq. 3.20
solveFlongODE(k,C,dt,index);
}
}
//********* total force *********//
// sum and convert to N
for (int i=0; i<3; i++) // loop over dimensions
{
ParmarBassetForce[i] = Fshort[i];
for (int k=0; k<2; k++) // loop over F1, F2
ParmarBassetForce[i] += FHist_[k][0][index][i];
}
ParmarBassetForce *= newton;
// Set value fields and write the probe
if(probeIt_)
{
scalar ReRef = 0.75/(gH0_[index][0]-0.105);
vector Flong1(0,0,0);
vector Flong2(0,0,0);
for (int i=0; i<3; i++) // loop over dimensions
{
Flong1[i] = FHist_[0][0][index][i];
Flong2[i] = FHist_[1][0][index][i];
}
//
// relative velocity (m/s)
vector UrelNoSmooth;
vector ddtUrelNoSmooth;
if(forceSubM(0).interpolation())
UrelNoSmooth = UrelNoSmoothInterpolator_.interpolate(position,cellI);
else
UrelNoSmooth = UrelNoSmooth_[cellI];
// acceleration (m/s2)
if(forceSubM(0).interpolation())
ddtUrelNoSmooth = ddtUrelNoSmoothInterpolator_.interpolate(position,cellI);
else
ddtUrelNoSmooth = ddtUrelNoSmooth_[cellI];
UrelNoSmooth /= mps;
ddtUrelNoSmooth /= mpss;
//
#include "setupProbeModelfields.H"
vValues.append(ParmarBassetForce); //first entry must the be the force
vValues.append(Urel);
vValues.append(ddtUrel);
//
vValues.append(UrelNoSmooth);
vValues.append(ddtUrelNoSmooth);
//
vValues.append(Fshort);
vValues.append(Flong1);
vValues.append(Flong2);
sValues.append(ReRef);
sValues.append(tRef_[index][0]);
sValues.append(mRef_[index][0]);
sValues.append(lRef_[index][0]);
sValues.append(r);
sValues.append(t0);
sValues.append(nKnown);
particleCloud_.probeM().writeProbe(index, sValues, vValues);
}
}
else
{
// not on CPU
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
// write particle based data to global array
forceSubM(0).partToArray(index,ParmarBassetForce,vector::zero);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::reAllocArrays() const
{
Pout << "ParmarBassetForce::reAllocArrays..." << endl;
particleCloud_.dataExchangeM().allocateArray(gH0_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(tRef_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(mRef_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(lRef_, NOTONCPU,1);
for (int i=0; i<nHist_; i++)
{
particleCloud_.dataExchangeM().allocateArray(ddtUrelHist_[i],NOTONCPU,3);
particleCloud_.dataExchangeM().allocateArray(rHist_ [i],NOTONCPU,1);
}
for (int i=0; i<2*discOrder_+1; i++)
for (int k=0; k<2; k++)
particleCloud_.dataExchangeM().allocateArray(FHist_[k][i],0.0,3);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
{
scalar cbrtr = cbrt(r); // cube root of r
scalar gamma = cbrtr*pow(dt,0.25);
/*
scalar K0 = 2./(9.*pow(r,0.666)) *
(
6.*gamma*gamma/(gamma*gamma*gamma+1)
+ 2.*sqrt(3.)*atan(sqrt(3.)*gamma/(2.-gamma))
+ log((gamma*gamma-gamma+1.)/(gamma+1.)/(gamma+1.))
); // int_0^dt K(r,xi) dxi
*/
scalar K0 = 2./(9.*cbrtr*cbrtr) *
(
- 0.37224 * gamma
+ 12.1587 * gamma*gamma
- 6.4884 * gamma*gamma*gamma
); // third order approximation to the eq. above
return K0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar Foam::ParmarBassetForce::trapWeight(int i, int n) const
{
if ( (i==1) || (i==(n-1)) )
return 0.5;
else
return 1.0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_ddtUrelHist(const vector& ddtUrel, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
for (int j=nHist_-1; j>0; j--) // loop over past times
ddtUrelHist_[j][index][i] = ddtUrelHist_[j-1][index][i];
ddtUrelHist_[0][index][i] = ddtUrel[i];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_rHist(scalar r, int index) const
{
for (int j=nHist_-1; j>0; j--) // loop over past times
rHist_[j][index][0] = rHist_[j-1][index][0];
rHist_[0][index][0] = r;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_FHist(const vector& F1, const vector& F2, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=FHistSize_-1; j>0; j--) // loop over past times
FHist_[k][j][index][i] = FHist_[k][j-1][index][i];
FHist_[0][0][index][i] = F1[i];
FHist_[1][0][index][i] = F2[i];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const
{
for (int j=0; j<4; j++) // loop over terms
{
C[j] = (c[k][j][0] * pow((t0 + c[k][j][1]),c[k][j][2])); //t0 dependency
C[j] *= (1 + chi[k][j][0]*(r-1) + chi[k][j][1]*(r-1)*(r-1)); // r dependency
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int index) const
{
if (discOrder_==1)
{
for (int i=0; i<3; i++) // loop over dimensions
{
FHist_[k][0][index][i] =
(
( C[1]/dt+2/(dt*dt)) * FHist_[k][1][index][i]
- ( 1/(dt*dt)) * FHist_[k][2][index][i]
- (C[2]+C[3]/dt ) * ddtUrelHist_[nInt_ ][index][i]
+ ( C[3]/dt ) * ddtUrelHist_[nInt_+1][index][i]
) / (C[0]+C[1]/dt+1/(dt*dt)); // Eq. 3.20 using first order temporal discretisation
}
}
if (discOrder_==2)
{
for (int i=0; i<3; i++) // loop over dimensions
{
FHist_[k][0][index][i] =
(
( 4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[k][1][index][i]
- ( C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[k][2][index][i]
+ ( 8/(4*dt*dt)) * FHist_[k][3][index][i]
- ( 1/(4*dt*dt)) * FHist_[k][4][index][i]
- (C[2] + 3*C[3]/(2*dt) ) * ddtUrelHist_[nInt_ ][index][i]
+ ( 4*C[3]/(2*dt) ) * ddtUrelHist_[nInt_+1][index][i]
- ( C[3]/(2*dt) ) * ddtUrelHist_[nInt_+2][index][i]
) / (C[0] + 3*C[1]/(2*dt) + 9/(4*dt*dt)); // Eq. 3.20 using second order temporal discretisation
}
}
}
void Foam::ParmarBassetForce::rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
// rescale ddtU history
for (int j=0; j<nHist_; j++) // loop over past times
if (ddtUrelHist_[j][index][i] != NOTONCPU)
ddtUrelHist_[j][index][i] *= lScale/(tScale*tScale);
// rescale F1, F2 history
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=0; j<FHistSize_; j++) // loop over past times
FHist_[k][j][index][i] *= mScale*lScale/(tScale*tScale);
}
// rescale r history
for (int j=0; j<nHist_; j++) // loop over past times
if (rHist_[j][index][0] != NOTONCPU)
rHist_[j][index][0] /= pow(rScale,1.5);
}
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
Class
ParmarBassetForce
SourceFiles
ParmarBassetForce.C
\*---------------------------------------------------------------------------*/
#ifndef ParmarBassetForce_H
#define ParmarBassetForce_H
#include "forceModel.H"
#include "dataExchangeModel.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ParmarBassetForce Declaration
\*---------------------------------------------------------------------------*/
class ParmarBassetForce
:
public forceModel
{
private:
dictionary propsDict_;
word velFieldName_;
const volVectorField& U_;
word UsFieldName_;
const volVectorField& Us_;
label nInt_; //no of timesteps to solve full integral
label discOrder_; //ODE discretisation order
label nHist_; //no of timesteps to save ddtUrel history for
label FHistSize_;
mutable List<double**> ddtUrelHist_;
mutable List<double**> rHist_;
mutable List<List<double**>> FHist_;
mutable double** gH0_;
mutable double** tRef_;
mutable double** mRef_;
mutable double** lRef_;
mutable volVectorField Urel_;
mutable volVectorField ddtUrel_;
autoPtr<smoothingModel> smoothingModel_;
private:
scalar calculateK0(scalar r, scalar dt) const;
scalar trapWeight(int i, int n) const;
void update_ddtUrelHist(const vector& ddtUrel, int index) const;
void update_rHist(scalar r, int index) const;
void update_FHist(const vector& F1, const vector& F2, int index) const;
void calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const;
void solveFlongODE(int k, double C[4], scalar dt, int index) const;
void rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const;
public:
//- Runtime type information
TypeName("ParmarBassetForce");
// Constructors
//- Construct from components
ParmarBassetForce
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~ParmarBassetForce();
// Member Functions
void setForce() const;
void reAllocArrays() const;
inline const smoothingModel& smoothingM() const
{
return smoothingModel_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -97,14 +97,16 @@ gradPForceSmooth::gradPForceSmooth
if (modelType_ == "B") if (modelType_ == "B")
{ {
FatalError <<"using model gradPForceSmooth with model type B is not valid\n" << abort(FatalError); FatalError <<"using model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
}else if (modelType_ == "Bfull") }
else if (modelType_ == "Bfull")
{ {
if(forceSubM(0).switches()[1]) if(forceSubM(0).switches()[1])
{ {
Info << "Using treatForceDEM false!" << endl; Info << "Using treatForceDEM false!" << endl;
forceSubM(0).setSwitches(1,false); // treatForceDEM = false forceSubM(0).setSwitches(1,false); // treatForceDEM = false
} }
}else // modelType_=="A" }
else // modelType_=="A"
{ {
if (!forceSubM(0).switches()[1]) if (!forceSubM(0).switches()[1])
{ {
@ -113,7 +115,9 @@ gradPForceSmooth::gradPForceSmooth
} }
} }
if (propsDict_.found("useU")) useU_=true; if (propsDict_.found("useU"))
useU_=true;
if (propsDict_.found("useAddedMass")) if (propsDict_.found("useAddedMass"))
{ {
addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass")); addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass"));
@ -131,15 +135,14 @@ gradPForceSmooth::gradPForceSmooth
particleCloud_.probeM().scalarFields_.append("Vs"); particleCloud_.probeM().scalarFields_.append("Vs");
particleCloud_.probeM().scalarFields_.append("rho"); particleCloud_.probeM().scalarFields_.append("rho");
particleCloud_.probeM().writeHeader(); particleCloud_.probeM().writeHeader();
pSmooth_ = p_;
} }
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
gradPForceSmooth::~gradPForceSmooth() gradPForceSmooth::~gradPForceSmooth()
{} {
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -159,23 +162,14 @@ void gradPForceSmooth::setForce() const
volScalarField pFull = pSmooth_ + rho_*gh_; volScalarField pFull = pSmooth_ + rho_*gh_;
gradPField = fvc::grad(pFull); gradPField = fvc::grad(pFull);
}
}else{ else
{
smoothingM().smoothen(pSmooth_); smoothingM().smoothen(pSmooth_);
gradPField = fvc::grad(pSmooth_); gradPField = fvc::grad(pSmooth_);
} }
/*if (useU_)
{
// const volScalarField& voidfraction_ = particleCloud_.mesh().lookupObject<volScalarField> ("voidfraction");
volScalarField U2 = U_&U_;// *voidfraction_*voidfraction_;
if (useRho_)
gradPField = fvc::grad(0.5*U2);
else
gradPField = fvc::grad(0.5*forceSubM(0).rhoField()*U2);
}*/
vector gradP; vector gradP;
scalar Vs; scalar Vs;
scalar rho; scalar rho;
@ -200,7 +194,8 @@ void gradPForceSmooth::setForce() const
if (forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!) if (forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!)
{ {
gradP = gradPInterpolator_.interpolate(position,cellI); gradP = gradPInterpolator_.interpolate(position,cellI);
}else }
else
{ {
gradP = gradPField[cellI]; gradP = gradPField[cellI];
} }

View File

@ -67,11 +67,37 @@ virtualMassForce::virtualMassForce
propsDict_(dict.subDict(typeName + "Props")), propsDict_(dict.subDict(typeName + "Props")),
velFieldName_(propsDict_.lookup("velFieldName")), velFieldName_(propsDict_.lookup("velFieldName")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)), U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
UsFieldName_(propsDict_.lookup("granVelFieldName")),
Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
phiFieldName_(propsDict_.lookup("phiFieldName")), phiFieldName_(propsDict_.lookup("phiFieldName")),
phi_(sm.mesh().lookupObject<surfaceScalarField> (phiFieldName_)), phi_(sm.mesh().lookupObject<surfaceScalarField> (phiFieldName_)),
UrelOldRegName_(typeName + "UrelOld"), UrelOldRegName_(typeName + "UrelOld"),
splitUrelCalculation_(propsDict_.lookupOrDefault<bool>("splitUrelCalculation",false)), splitUrelCalculation_(propsDict_.lookupOrDefault<bool>("splitUrelCalculation",false)),
Cadd_(0.5) useUs_(propsDict_.lookupOrDefault<bool>("useUs",false)),
useFelderhof_(propsDict_.lookupOrDefault<bool>("useFelderhof",false)),
Cadd_(0.5),
DDtUrel_
( IOobject
(
"DDtUrel",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(0,1,-2,0,0,0,0), vector(0,0,0))
),
smoothingModel_
(
smoothingModel::New
(
propsDict_,
sm
)
)
{ {
particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false); particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false);
@ -93,17 +119,36 @@ virtualMassForce::virtualMassForce
Cadd_ = readScalar(propsDict_.lookup("Cadd")); Cadd_ = readScalar(propsDict_.lookup("Cadd"));
Info << "Virtual mass model: using non-standard Cadd = " << Cadd_ << endl; Info << "Virtual mass model: using non-standard Cadd = " << Cadd_ << endl;
} }
if(useUs_)
{
Info << "Virtual mass model: using averaged Us \n";
Info << "WARNING: ignoring virtual mass of relative particle motion/collisions \n";
if(splitUrelCalculation_)
{
FatalError << "Virtual mass model: useUs=true requires splitUrelCalculation_=false" << abort(FatalError);
}
}
if(useFelderhof_)
{
Info << "Virtual mass model: using Cadd correlation by Felderhof \n";
Info << "WARNING: ignoring user-set Cadd \n";
if(!dict.lookupOrDefault<bool>("getParticleDensities",false))
FatalError << "Virtual mass model: useFelderhof=true requires getParticleDensities=true in couplingProperties" << abort(FatalError);
}
particleCloud_.checkCG(true); particleCloud_.checkCG(true);
//Append the field names to be probed //Append the field names to be probed
particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().initialize(typeName, typeName+".logDat");
particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force
particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().vectorFields_.append("acceleration");
particleCloud_.probeM().vectorFields_.append("UrelOld"); particleCloud_.probeM().scalarFields_.append("Cadd");
particleCloud_.probeM().vectorFields_.append("ddtUrel");
particleCloud_.probeM().scalarFields_.append("Vs"); particleCloud_.probeM().scalarFields_.append("Vs");
particleCloud_.probeM().scalarFields_.append("rho"); particleCloud_.probeM().scalarFields_.append("rho");
particleCloud_.probeM().scalarFields_.append("voidfraction");
particleCloud_.probeM().scalarFields_.append("specificGravity");
particleCloud_.probeM().writeHeader(); particleCloud_.probeM().writeHeader();
} }
@ -122,18 +167,18 @@ void virtualMassForce::setForce() const
scalar dt = U_.mesh().time().deltaT().value(); scalar dt = U_.mesh().time().deltaT().value();
vector position(0,0,0); //Compute acceleration field
vector Ufluid(0,0,0);
vector Ur(0,0,0);
vector DDtU(0,0,0);
//Compute extra vfields in case it is needed
volVectorField DDtU_(0.0*U_/U_.mesh().time().deltaT());
if(splitUrelCalculation_) if(splitUrelCalculation_)
DDtU_ = fvc::ddt(U_) + fvc::div(phi_, U_); //Total Derivative of fluid velocity DDtUrel_ = fvc::ddt(U_) + fvc::div(phi_, U_); //Total Derivative of fluid velocity
else if(useUs_)
DDtUrel_ = fvc::ddt(U_) + fvc::div(phi_, U_) - fvc::ddt(Us_); //Total Derivative of fluid velocity minus average particle velocity
// smoothen
smoothingM().smoothen(DDtUrel_);
interpolationCellPoint<vector> UInterpolator_(U_); interpolationCellPoint<vector> UInterpolator_(U_);
interpolationCellPoint<vector> DDtUInterpolator_(DDtU_); interpolationCellPoint<vector> DDtUrelInterpolator_(DDtUrel_);
interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
#include "setupProbeModel.H" #include "setupProbeModel.H"
@ -142,37 +187,42 @@ void virtualMassForce::setForce() const
for(int index = 0;index < particleCloud_.numberOfParticles(); index++) for(int index = 0;index < particleCloud_.numberOfParticles(); index++)
{ {
vector virtualMassForce(0,0,0); vector virtualMassForce(0,0,0);
vector position(0,0,0);
vector DDtUrel(0,0,0);
scalar voidfraction(1);
scalar sg(1);
label cellI = particleCloud_.cellIDs()[index][0]; label cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1) // particle Found if (cellI > -1) // particle Found
{ {
// particle position
if(forceSubM(0).interpolation()) if(forceSubM(0).interpolation())
{
position = particleCloud_.position(index); position = particleCloud_.position(index);
Ufluid = UInterpolator_.interpolate(position,cellI);
} //********* acceleration value *********//
else if(splitUrelCalculation_ || useUs_ )
{ {
Ufluid = U_[cellI]; // DDtUrel from acceleration field
}
if(splitUrelCalculation_) //if split, just use total derivative of fluid velocity
if(forceSubM(0).interpolation()) if(forceSubM(0).interpolation())
{ DDtUrel = DDtUrelInterpolator_.interpolate(position,cellI);
DDtU = DDtUInterpolator_.interpolate(position,cellI); else
DDtUrel = DDtUrel_[cellI];
} }
else else
{ {
DDtU = DDtU_[cellI]; // DDtUrel from UrelOld
} vector Ufluid(0,0,0);
else
{
vector Us = particleCloud_.velocity(index);
Ur = Ufluid - Us;
}
// relative velocity
if(forceSubM(0).interpolation())
Ufluid = UInterpolator_.interpolate(position,cellI);
else
Ufluid = U_[cellI];
vector Us = particleCloud_.velocity(index);
vector Urel = Ufluid - Us;
//Check of particle was on this CPU the last step //Check of particle was on this CPU the last step
if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step
@ -185,19 +235,46 @@ void virtualMassForce::setForce() const
for(int j=0;j<3;j++) for(int j=0;j<3;j++)
{ {
UrelOld[j] = UrelOld_[index][j]; UrelOld[j] = UrelOld_[index][j];
UrelOld_[index][j] = Ur[j]; UrelOld_[index][j] = Urel[j];
} }
if(haveUrelOld_ ) //only compute force if we have old (relative) velocity if(haveUrelOld_ ) //only compute force if we have old (relative) velocity
ddtUrel = (Ur-UrelOld)/dt; ddtUrel = (Urel-UrelOld)/dt;
}
if(splitUrelCalculation_) //we can always compute the total derivative in case we split
ddtUrel = DDtU;
//********* Cadd value *********//
scalar rho = forceSubM(0).rhoField()[cellI];
scalar ds = 2*particleCloud_.radius(index); scalar ds = 2*particleCloud_.radius(index);
scalar Vs = ds*ds*ds*M_PI/6; scalar Vs = ds*ds*ds*M_PI/6;
scalar rho = forceSubM(0).rhoField()[cellI];
virtualMassForce = Cadd_ * rho * Vs * ddtUrel; scalar Cadd;
if (useFelderhof_)
{
scalar epsilons(0);
scalar logsg(0);
if(forceSubM(0).interpolation())
voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
else
voidfraction = voidfraction_[cellI];
sg = particleCloud_.particleDensity(index) / rho;
logsg = log(sg);
epsilons = 1-voidfraction;
Cadd = 0.5
+ ( 0.047*logsg + 0.13)*epsilons
+ (-0.066*logsg - 0.58)*epsilons*epsilons
+ ( 1.42)*epsilons*epsilons*epsilons;
}
else
{
// use predefined value
Cadd = Cadd_;
}
//********* calculate force *********//
virtualMassForce = Cadd_ * rho * Vs * DDtUrel;
if( forceSubM(0).verbose() ) //&& index>100 && index < 105) if( forceSubM(0).verbose() ) //&& index>100 && index < 105)
{ {
@ -210,15 +287,17 @@ void virtualMassForce::setForce() const
{ {
#include "setupProbeModelfields.H" #include "setupProbeModelfields.H"
vValues.append(virtualMassForce); //first entry must the be the force vValues.append(virtualMassForce); //first entry must the be the force
vValues.append(Ur); vValues.append(DDtUrel);
vValues.append(UrelOld); sValues.append(Cadd);
vValues.append(ddtUrel);
sValues.append(Vs); sValues.append(Vs);
sValues.append(rho); sValues.append(rho);
sValues.append(voidfraction);
sValues.append(sg);
particleCloud_.probeM().writeProbe(index, sValues, vValues); particleCloud_.probeM().writeProbe(index, sValues, vValues);
} }
} }
else //particle not on this CPU else //particle not on this CPU
if (!useUs_ || !splitUrelCalculation_ )
UrelOld_[index][0] = NOTONCPU; UrelOld_[index][0] = NOTONCPU;
// write particle based data to global array // write particle based data to global array

View File

@ -67,6 +67,14 @@ private:
const volVectorField& U_; const volVectorField& U_;
word voidfractionFieldName_;
const volScalarField& voidfraction_;
word UsFieldName_;
const volVectorField& Us_;
word phiFieldName_; word phiFieldName_;
const surfaceScalarField& phi_; const surfaceScalarField& phi_;
@ -76,7 +84,15 @@ private:
const bool splitUrelCalculation_; //indicator to split calculation of Urel between CFDEM and LIGGGHTS const bool splitUrelCalculation_; //indicator to split calculation of Urel between CFDEM and LIGGGHTS
//requires the integration fix to take dv/dt into account! //requires the integration fix to take dv/dt into account!
mutable double Cadd_; //virtual mass constant mutable bool useUs_; //indicator to use averaged Us instead of v_p
mutable bool useFelderhof_; //indicator to use Cadd values from Felderhof
mutable double Cadd_; //custom Cadd value
mutable volVectorField DDtUrel_;
autoPtr<smoothingModel> smoothingModel_;
public: public:
@ -100,6 +116,11 @@ public:
// Member Functions // Member Functions
void setForce() const; void setForce() const;
inline const smoothingModel& smoothingM() const
{
return smoothingModel_;
}
}; };

View File

@ -147,7 +147,7 @@ void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict)
// if this makes troubles try floor((startTime_+SMALL)/.. as above // if this makes troubles try floor((startTime_+SMALL)/.. as above
firstCouplingStep_ = floor((startTime_+SMALL-simStartTime)/DEMts/couplingInterval); firstCouplingStep_ = floor((startTime_+SMALL-simStartTime)/DEMts/couplingInterval);
lastCouplingStep_ = floor((endTime_+SMALL-simStartTime)/DEMts/couplingInterval); lastCouplingStep_ = floor((endTime_+SMALL-simStartTime)/DEMts/couplingInterval);
couplingStepInterval_ = floor(timeInterval_+SMALL/DEMts/couplingInterval); couplingStepInterval_ = floor((timeInterval_+SMALL)/DEMts/couplingInterval);
} }
else //runEveryCouplingStep or writeStep else //runEveryCouplingStep or writeStep
{ {
@ -279,6 +279,11 @@ void liggghtsCommandModel::parseCommandList(wordList& commandList,labelList& lab
else if (add=="dot") add = "."; else if (add=="dot") add = ".";
else if (add=="dotdot") add = ".."; else if (add=="dotdot") add = "..";
else if (add=="slash") add = "/"; else if (add=="slash") add = "/";
else if (add=="dollar") add = "$";
else if (add=="curlyOpen") add = "{";
else if (add=="curlyClose") add = "}";
else if (add=="squareOpen") add = "[";
else if (add=="squareClose") add = "]";
else if (add=="noBlanks") // no blanks after the following words else if (add=="noBlanks") // no blanks after the following words
{ {
add = ""; add = "";

View File

@ -0,0 +1,484 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "massTransferGunn.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(massTransferGunn, 0);
addToRunTimeSelectionTable(massTransferModel, massTransferGunn, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
massTransferGunn::massTransferGunn
(
const dictionary& dict,
cfdemCloudEnergy& sm
)
:
massTransferModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
multiTypes_(false),
expSherwood_(propsDict_.lookupOrDefault<bool>("expSherwood",false)),
interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
SPartFluidName_(propsDict_.lookupOrDefault<word>("SPartFluidName","SPartFluid")),
SPartFluid_
( IOobject
(
SPartFluidName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
),
SPartFluidCoeffName_(propsDict_.lookupOrDefault<word>("SPartFluidCoeffName","SPartFluidCoeff")),
SPartFluidCoeff_
( IOobject
(
SPartFluidCoeffName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,-1,0,0,0,0), 0.0)
),
ReField_
( IOobject
(
"ReField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
),
ShField_
( IOobject
(
"ShField",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
),
concFieldName_(propsDict_.lookupOrDefault<word>("concFieldName","C")),
concField_(sm.mesh().lookupObject<volScalarField> (concFieldName_)),
satConcFieldName_(propsDict_.lookupOrDefault<word>("satConcFieldName","Cs")),
satConcField_(sm.mesh().lookupObject<volScalarField> (satConcFieldName_)),
voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
maxSource_(1e30),
velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
coupleDEM_(propsDict_.lookupOrDefault<bool>("coupleDEM",false)),
partMassFluxName_(propsDict_.lookupOrDefault<word>("partMassFluxName","convectiveMassFlux")),
partMassFluxCoeffRegName_(typeName + "partMassFluxCoeff"),
partReRegName_(typeName + "partRe"),
partShRegName_(typeName + "partSh"),
scaleDia_(1.),
typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0)))
{
particleCloud_.registerParticleProperty<double**>(partMassFluxName_,1);
if(implicit_)
{
particleCloud_.registerParticleProperty<double**>(partMassFluxCoeffRegName_,1);
}
if(verbose_)
{
particleCloud_.registerParticleProperty<double**>(partReRegName_,1);
particleCloud_.registerParticleProperty<double**>(partShRegName_,1);
}
if (propsDict_.found("maxSource"))
{
maxSource_=readScalar(propsDict_.lookup ("maxSource"));
Info << "limiting eulerian source field to: " << maxSource_ << endl;
}
if (!implicit_)
{
SPartFluidCoeff_.writeOpt() = IOobject::NO_WRITE;
}
if (verbose_)
{
ReField_.writeOpt() = IOobject::AUTO_WRITE;
ShField_.writeOpt() = IOobject::AUTO_WRITE;
ReField_.write();
ShField_.write();
if (expSherwood_)
{
FatalError <<"Cannot read and create ShField at the same time!\n" << abort(FatalError);
}
}
if (propsDict_.found("scale") && typeCG_.size()==1)
{
// if "scale" is specified and there's only one single type, use "scale"
scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
typeCG_[0] = scaleDia_;
}
else if (typeCG_.size()>1)
{
multiTypes_ = true;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
massTransferGunn::~massTransferGunn()
{
}
// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
void massTransferGunn::calcMassContribution()
{
double**& partMassFlux_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxName_);
// reset Scalar field
SPartFluid_.primitiveFieldRef() = 0.0;
#ifdef compre
const volScalarField mufField = particleCloud_.turbulence().mu();
#else
const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
#endif
const volScalarField& D0Field_ = D0Field();
if (typeCG_.size()>1 || typeCG_[0] > 1)
{
Info << "massTransferGunn using scale = " << typeCG_ << endl;
}
else if (particleCloud_.cg() > 1)
{
scaleDia_=particleCloud_.cg();
Info << "massTransferGunn using scale from liggghts cg = " << scaleDia_ << endl;
}
// calc La based heat flux
scalar voidfraction(1);
vector Ufluid(0,0,0);
scalar Cfluid(0);
scalar Csfluid(0);
label cellI=0;
vector Us(0,0,0);
scalar ds(0);
scalar ds_scaled(0);
scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
scalar muf(0);
scalar rhof(0);
scalar magUr(0);
scalar Rep(0);
scalar Sc(0);
scalar Shp(0);
scalar cg = scaleDia_;
label partType = 1;
interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
interpolationCellPoint<vector> UInterpolator_(U_);
interpolationCellPoint<scalar> CInterpolator_(concField_);
for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
{
cellI = particleCloud_.cellIDs()[index][0];
if(cellI >= 0)
{
if(interpolation_)
{
vector position = particleCloud_.position(index);
voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
Ufluid = UInterpolator_.interpolate(position,cellI);
Cfluid = CInterpolator_.interpolate(position,cellI);
}
else
{
voidfraction = voidfraction_[cellI];
Ufluid = U_[cellI];
Cfluid = concField_[cellI];
}
if (voidfraction < 0.01)
voidfraction = 0.01;
if (multiTypes_)
{
partType = particleCloud_.particleType(index);
cg = typeCG_[partType - 1];
scaleDia3 = cg*cg*cg;
}
// calc relative velocity
Us = particleCloud_.velocity(index);
magUr = mag(Ufluid - Us);
ds = 2.*particleCloud_.radius(index);
ds_scaled = ds/cg;
muf = mufField[cellI];
rhof = rho_[cellI];
Csfluid = satConcField_[cellI];
Rep = ds_scaled * magUr * voidfraction * rhof/ muf;
scalar D0 = D0Field_[cellI];
Sc = max(SMALL, muf / (rhof*D0));
Shp = Sherwood(voidfraction, Rep, Sc);
scalar h = D0 * Shp / ds_scaled;
scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere
// calc convective heat flux [W]
massFlux(index, h, As, Cfluid, Csfluid, scaleDia3);
if(verbose_)
{
double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
double**& partSh_ = particleCloud_.getParticlePropertyRef<double**>(partShRegName_);
partRe_[index][0] = Rep;
partSh_[index][0] = Shp;
}
if(verbose_ && index >=0 && index <2)
{
Pout << "partMassFlux = " << partMassFlux_[index][0] << endl;
Pout << "magUr = " << magUr << endl;
Pout << "D0 = " << D0 << endl;
Pout << "rho = " << rhof << endl;
Pout << "h = " << h << endl;
Pout << "ds = " << ds << endl;
Pout << "ds_scaled = " << ds_scaled << endl;
Pout << "As = " << As << endl;
Pout << "muf = " << muf << endl;
Pout << "Rep = " << Rep << endl;
Pout << "Sc = " << Sc << endl;
Pout << "Shp = " << Shp << endl;
Pout << "voidfraction = " << voidfraction << endl;
Pout << "Cfluid = " << Cfluid << endl;
Pout << "Csfluid = " << Csfluid << endl;
}
}
}
particleCloud_.averagingM().setScalarSum
(
SPartFluid_,
partMassFlux_,
particleCloud_.particleWeights(),
NULL
);
SPartFluid_.primitiveFieldRef() /= -SPartFluid_.mesh().V();
if(implicit_)
{
double**& partMassFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxCoeffRegName_);
SPartFluidCoeff_.primitiveFieldRef() = 0.0;
particleCloud_.averagingM().setScalarSum
(
SPartFluidCoeff_,
partMassFluxCoeff_,
particleCloud_.particleWeights(),
NULL
);
SPartFluidCoeff_.primitiveFieldRef() /= -SPartFluidCoeff_.mesh().V();
}
if(verbose_)
{
double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
double**& partSh_ = particleCloud_.getParticlePropertyRef<double**>(partShRegName_);
ReField_.primitiveFieldRef() = 0.0;
ShField_.primitiveFieldRef() = 0.0;
particleCloud_.averagingM().resetWeightFields();
particleCloud_.averagingM().setScalarAverage
(
ReField_,
partRe_,
particleCloud_.particleWeights(),
particleCloud_.averagingM().UsWeightField(),
NULL
);
particleCloud_.averagingM().resetWeightFields();
particleCloud_.averagingM().setScalarAverage
(
ShField_,
partSh_,
particleCloud_.particleWeights(),
particleCloud_.averagingM().UsWeightField(),
NULL
);
}
// limit source term in explicit treatment
if(!implicit_)
{
forAll(SPartFluid_,cellI)
{
scalar EuFieldInCell = SPartFluid_[cellI];
if(mag(EuFieldInCell) > maxSource_ )
{
Pout << "limiting source term\n" << endl ;
SPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
}
}
}
SPartFluid_.correctBoundaryConditions();
volScalarField minParticleWeights = particleCloud_.averagingM().UsWeightField();
Info << "Minimum Particle Weight " << gMin(minParticleWeights) << endl;
Info << "Minimum Fluid Concentration: " << gMin(concField_) << endl;
Info << "Maximum Fluid Concentration: " << gMax(concField_) << endl;
}
void massTransferGunn::addMassContribution(volScalarField& Sm) const
{
Sm += SPartFluid_;
}
void massTransferGunn::addMassCoefficient(volScalarField& Smi) const
{
if(implicit_)
{
Smi += SPartFluidCoeff_;
}
}
scalar massTransferGunn::Sherwood(scalar voidfraction, scalar Rep, scalar Sc) const
{
return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
(1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Sc,0.33)) +
(1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) *
Foam::pow(Rep,0.7) * Foam::pow(Sc,0.33);
}
void massTransferGunn::massFlux(label index, scalar h, scalar As, scalar Cfluid, scalar Csfluid, scalar cg3)
{
scalar hAs = h * As * cg3;
double**& partMassFlux_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxName_);
if (particleCloud_.getParticleEffVolFactors())
{
scalar effVolFac = particleCloud_.particleEffVolFactor(index);
hAs *= effVolFac;
}
partMassFlux_[index][0] = - hAs * Csfluid;
if(!implicit_)
{
partMassFlux_[index][0] += hAs * Cfluid;
}
else
{
double**& partMassFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxCoeffRegName_);
partMassFluxCoeff_[index][0] = hAs;
}
}
void massTransferGunn::giveData()
{
Info << "total convective particle-fluid mass flux [kg/s] (Eulerian) = " << gSum(SPartFluid_*1.0*SPartFluid_.mesh().V()) << endl;
double**& partMassFlux_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxName_);
particleCloud_.dataExchangeM().giveData(partMassFluxName_,"scalar-atom", partMassFlux_);
}
void massTransferGunn::postFlow()
{
// send mass flux to DEM
if (coupleDEM_)
{
if(implicit_)
{
label cellI;
scalar Cfluid(0.0);
interpolationCellPoint<scalar> CInterpolator_(concField_);
double**& partMassFlux_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxName_);
double**& partMassFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partMassFluxCoeffRegName_);
for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
{
cellI = particleCloud_.cellIDs()[index][0];
if(cellI >= 0)
{
if(interpolation_)
{
vector position = particleCloud_.position(index);
Cfluid = CInterpolator_.interpolate(position,cellI);
}
else
{
Cfluid = concField_[cellI];
}
partMassFlux_[index][0] += Cfluid * partMassFluxCoeff_[index][0];
if(verbose_ && index >=0 && index <2)
{
Pout << "partMassFlux = " << partMassFlux_[index][0] << endl;
Pout << "Cfluid = " << Cfluid << endl;
}
}
}
}
giveData();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
Correlation for Sherwood number according to
Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
\*---------------------------------------------------------------------------*/
#ifndef massTransferGunn_H
#define massTransferGunn_H
#include "fvCFD.H"
#include "cfdemCloudEnergy.H"
#include "massTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class massTransferGunn Declaration
\*---------------------------------------------------------------------------*/
class massTransferGunn
:
public massTransferModel
{
protected:
dictionary propsDict_;
bool multiTypes_;
bool expSherwood_;
bool interpolation_;
bool verbose_;
bool implicit_;
word SPartFluidName_;
volScalarField SPartFluid_;
word SPartFluidCoeffName_;
volScalarField SPartFluidCoeff_;
volScalarField ReField_;
volScalarField ShField_;
word concFieldName_;
const volScalarField& concField_; // ref to concentration field
word satConcFieldName_;
const volScalarField& satConcField_; // ref to saturation concentration field
word voidfractionFieldName_;
const volScalarField& voidfraction_; // ref to voidfraction field
scalar maxSource_; // max (limited) value of src field
word velFieldName_;
const volVectorField& U_;
word densityFieldName_;
const volScalarField& rho_;
bool coupleDEM_;
const word partMassFluxName_;
const word partMassFluxCoeffRegName_;
const word partReRegName_;
const word partShRegName_;
mutable scalar scaleDia_;
scalarList typeCG_;
scalar Sherwood(scalar, scalar, scalar) const;
virtual void giveData();
virtual void massFlux(label, scalar, scalar, scalar, scalar, scalar cg3 = 1.0);
public:
//- Runtime type information
TypeName("massTransferGunn");
// Constructors
//- Construct from components
massTransferGunn
(
const dictionary& dict,
cfdemCloudEnergy& sm
);
// Destructor
virtual ~massTransferGunn();
// Member Functions
void addMassContribution(volScalarField&) const;
void addMassCoefficient(volScalarField&) const;
void calcMassContribution();
void postFlow();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "massTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(massTransferModel, 0);
defineRunTimeSelectionTable(massTransferModel, dictionary);
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
massTransferModel::massTransferModel
(
const dictionary& dict,
cfdemCloudEnergy& sm
)
:
dict_(dict),
particleCloud_(sm),
transportProperties_
(
IOobject
(
"transportProperties",
sm.mesh().time().constant(),
sm.mesh(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
D0Field_
(
IOobject
(
"D0",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.)
),
CsField_
(
IOobject
(
"Cs",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-3,0,0,0,0,0), 0.)
)
{
// build constant fields for single phase case
if (!particleCloud_.multiphase())
{
D0Field_ = volScalarField
(
IOobject
(
"D0",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("D"))
);
CsField_ = volScalarField
(
IOobject
(
"Cs",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("Cs"))
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
massTransferModel::~massTransferModel()
{}
// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * //
const volScalarField& massTransferModel::D0Field() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("D");
}
else
{
return D0Field_;
}
}
const volScalarField& massTransferModel::CsField() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("Cs");
}
else
{
return CsField_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,124 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#ifndef massTransferModel_H
#define massTransferModel_H
#include "fvCFD.H"
#include "cfdemCloudEnergy.H"
#include "probeModel.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class massTransferModel Declaration
\*---------------------------------------------------------------------------*/
class massTransferModel
{
protected:
// Protected data
const dictionary& dict_;
cfdemCloudEnergy& particleCloud_;
IOdictionary transportProperties_;
mutable volScalarField D0Field_;
mutable volScalarField CsField_;
public:
//- Runtime type information
TypeName("massTransferModel");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
massTransferModel,
dictionary,
(
const dictionary& dict,
cfdemCloudEnergy& sm
),
(dict,sm)
);
// Constructors
//- Construct from components
massTransferModel
(
const dictionary& dict,
cfdemCloudEnergy& sm
);
// Destructor
virtual ~massTransferModel();
// Selector
static autoPtr<massTransferModel> New
(
const dictionary& dict,
cfdemCloudEnergy& sm,
word massTransferType
);
// Member Functions
virtual void addMassContribution(volScalarField&) const = 0;
virtual void addMassCoefficient(volScalarField&) const = 0;
virtual void calcMassContribution() = 0;
virtual void postFlow() {}
virtual void solve() {}
const volScalarField& D0Field() const;
const volScalarField& CsField() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,67 @@
/*---------------------------------------------------------------------------*\
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 <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "massTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
autoPtr<massTransferModel> massTransferModel::New
(
const dictionary& dict,
cfdemCloudEnergy& sm,
word massTransferType
)
{
Info<< "Selecting massTransferModel "
<< massTransferType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(massTransferType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalError
<< "massTransferModel::New(const dictionary&, const spray&) : "
<< endl
<< " unknown massTransferModelType type "
<< massTransferType
<< ", constructor not in hash table" << endl << endl
<< " Valid massTransferModel types are :"
<< endl;
Info<< dictionaryConstructorTablePtr_->toc()
<< abort(FatalError);
}
return autoPtr<massTransferModel>(cstrIter()(dict,sm));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,265 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright (C) 2013- Graz University of
Technology, IPPT
-------------------------------------------------------------------------------
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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "constDiffAndTemporalSmoothing.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(constDiffAndTemporalSmoothing, 0);
addToRunTimeSelectionTable
(
smoothingModel,
constDiffAndTemporalSmoothing,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
constDiffAndTemporalSmoothing::constDiffAndTemporalSmoothing
(
const dictionary& dict,
cfdemCloud& sm
)
:
smoothingModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
lowerLimit_(readScalar(propsDict_.lookup("lowerLimit"))),
upperLimit_(readScalar(propsDict_.lookup("upperLimit"))),
smoothingLength_(dimensionedScalar("smoothingLength",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
gamma_(readScalar(propsDict_.lookup("smoothingStrength"))),
correctBoundary_(propsDict_.lookupOrDefault<bool>("correctBoundary",false)),
verbose_(false)
{
if(propsDict_.found("verbose"))
verbose_ = true;
if(propsDict_.found("smoothingLengthReferenceField"))
smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField")));
checkFields(sSmoothField_);
checkFields(vSmoothField_);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
constDiffAndTemporalSmoothing::~constDiffAndTemporalSmoothing()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool constDiffAndTemporalSmoothing::doSmoothing() const
{
return true;
}
void constDiffAndTemporalSmoothing::smoothen(volScalarField& fieldSrc) const
{
// Create scalar smooth field from virgin scalar smooth field template
volScalarField sSmoothField = sSmoothField_;
sSmoothField.dimensions().reset(fieldSrc.dimensions());
sSmoothField.ref()=fieldSrc.internalField();
sSmoothField.correctBoundaryConditions();
sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
sSmoothField.oldTime()=fieldSrc;
sSmoothField.oldTime().correctBoundaryConditions();
double deltaT = sSmoothField.mesh().time().deltaTValue();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
// do smoothing
solve
(
fvm::ddt(sSmoothField)
-fvm::laplacian(DT_, sSmoothField)
);
// bound sSmoothField_
forAll(sSmoothField,cellI)
{
sSmoothField[cellI]=max(lowerLimit_,min(upperLimit_,sSmoothField[cellI]));
}
// get data from working sSmoothField - will copy only values at new time
fieldSrc=sSmoothField;
fieldSrc.correctBoundaryConditions();
if(verbose_)
{
Info << "min/max(fieldoldTime) (unsmoothed): " << min(sSmoothField.oldTime()) << tab << max(sSmoothField.oldTime()) << endl;
Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void constDiffAndTemporalSmoothing::smoothen(volVectorField& fieldSrc) const
{
// Create scalar smooth field from virgin scalar smooth field template
volVectorField vSmoothField = vSmoothField_;
vSmoothField.dimensions().reset(fieldSrc.dimensions());
vSmoothField.ref()=fieldSrc.internalField();
vSmoothField.correctBoundaryConditions();
vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
vSmoothField.oldTime()=fieldSrc;
vSmoothField.oldTime().correctBoundaryConditions();
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT.value();
// do spacial smoothing
solve
(
fvm::ddt(vSmoothField)
-fvm::laplacian(DT_, vSmoothField)
);
// create fields for temporal smoothing
volVectorField refField = vSmoothField;
vSmoothField.ref()=fieldSrc.oldTime().internalField();
vSmoothField.correctBoundaryConditions();
vSmoothField.oldTime()=fieldSrc.oldTime();
vSmoothField.oldTime().correctBoundaryConditions();
// do temporal smoothing
solve
(
fvm::ddt(vSmoothField)
-
gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
);
// create temporally smoothened boundary field
if (correctBoundary_)
vSmoothField.boundaryFieldRef() = 1/(1+gamma_)*fieldSrc.oldTime().boundaryField() + gamma_/(1+gamma_)*fieldSrc.boundaryField();
// get data from working vSmoothField
fieldSrc=vSmoothField;
fieldSrc.correctBoundaryConditions();
if(verbose_)
{
Info << "min/max(fieldoldTime) (unsmoothed): " << min(vSmoothField.oldTime()) << tab << max(vSmoothField.oldTime()) << endl;
Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void constDiffAndTemporalSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
{
// Create scalar smooth field from virgin scalar smooth field template
volVectorField vSmoothField = vSmoothField_;
vSmoothField.dimensions().reset(fieldSrc.dimensions());
vSmoothField.ref()=fieldSrc.internalField();
vSmoothField.correctBoundaryConditions();
vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
vSmoothField.oldTime()=fieldSrc;
vSmoothField.oldTime().correctBoundaryConditions();
double sourceStrength = 1e5; //large number to keep reference values constant
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
DT_.value() = smoothingLengthReferenceField_.value()
* smoothingLengthReferenceField_.value() / deltaT.value();
tmp<volScalarField> NLarge
(
new volScalarField
(
IOobject
(
"NLarge",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
0.0
)
);
//loop over particles and map max particle diameter to Euler Grid
forAll(vSmoothField,cellI)
{
if ( mag(vSmoothField.oldTime().internalField()[cellI]) > 0.0f) // have a vector in the OLD vSmoothField, so keep it!
NLarge.ref()[cellI] = sourceStrength;
}
// do the smoothing
solve
(
fvm::ddt(vSmoothField)
-fvm::laplacian( DT_, vSmoothField)
==
NLarge() / deltaT * vSmoothField.oldTime() //add source to keep cell values constant
-fvm::Sp( NLarge() / deltaT, vSmoothField) //add sink to keep cell values constant
);
// get data from working vSmoothField
fieldSrc=vSmoothField;
fieldSrc.correctBoundaryConditions();
if(verbose_)
{
Info << "min/max(fieldoldTime) (unsmoothed): " << min(vSmoothField.oldTime()) << tab << max(vSmoothField.oldTime()) << endl;
Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,114 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright (C) 2013- Graz University of
Technology, IPPT
-------------------------------------------------------------------------------
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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
Class
constDiffAndTemporalSmoothing
SourceFiles
constDiffAndTemporalSmoothing.C
\*---------------------------------------------------------------------------*/
#ifndef constDiffAndTemporalSmoothing_H
#define constDiffAndTemporalSmoothing_H
#include "smoothingModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class constDiffAndTemporalSmoothing Declaration
\*---------------------------------------------------------------------------*/
class constDiffAndTemporalSmoothing
:
public smoothingModel
{
private:
dictionary propsDict_;
scalar lowerLimit_;
scalar upperLimit_;
dimensionedScalar smoothingLength_;
dimensionedScalar smoothingLengthReferenceField_;
mutable dimensionedScalar DT_;
scalar gamma_;
bool correctBoundary_;
bool verbose_;
public:
//- Runtime type information
TypeName("constDiffAndTemporalSmoothing");
// Constructors
//- Construct from components
constDiffAndTemporalSmoothing
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~constDiffAndTemporalSmoothing();
// Member Functions
bool doSmoothing() const;
//void dSmoothing(volScalarField&) const;
void smoothen(volScalarField&) const;
void smoothen(volVectorField&) const;
void smoothenReferenceField(volVectorField&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -64,11 +64,15 @@ SyamlalThermCond::~SyamlalThermCond()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void SyamlalThermCond::calcThermCond() void SyamlalThermCond::calcThermCond()
{ {
const volScalarField& kf0Field_ = kf0Field();
forAll(thermCondField_,cellI) forAll(thermCondField_,cellI)
{ {
if (1-voidfraction_[cellI] < SMALL) thermCondField_[cellI] = kf0_.value(); scalar kf0 = kf0Field_[cellI];
if (1-voidfraction_[cellI] < SMALL) thermCondField_[cellI] = kf0;
else if (voidfraction_[cellI] < SMALL) thermCondField_[cellI] = 0.0; else if (voidfraction_[cellI] < SMALL) thermCondField_[cellI] = 0.0;
else thermCondField_[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf0_.value(); else thermCondField_[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf0;
} }
thermCondField_.correctBoundaryConditions(); thermCondField_.correctBoundaryConditions();

View File

@ -83,6 +83,7 @@ ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond()
void ZehnerSchluenderThermCond::calcThermCond() void ZehnerSchluenderThermCond::calcThermCond()
{ {
const volScalarField& kf0Field_ = kf0Field();
calcPartKsField(); calcPartKsField();
scalar A = 0.0; scalar A = 0.0;
scalar B = 0.0; scalar B = 0.0;
@ -98,12 +99,13 @@ void ZehnerSchluenderThermCond::calcThermCond()
if(voidfraction > 1.0 - SMALL || partKsField_[cellI] < SMALL) partKsField_[cellI] = 0.0; if(voidfraction > 1.0 - SMALL || partKsField_[cellI] < SMALL) partKsField_[cellI] = 0.0;
else else
{ {
A = partKsField_[cellI]/kf0_.value(); scalar kf0 = kf0Field_[cellI];
A = partKsField_[cellI]/kf0;
B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11);
InvOnemBoA = 1.0/(1.0 - B/A); InvOnemBoA = 1.0/(1.0 - B/A);
C = (A - 1) * InvOnemBoA * InvOnemBoA * B/A * log(A/B) - (B - 1) * InvOnemBoA - 0.5 * (B + 1); C = (A - 1) * InvOnemBoA * InvOnemBoA * B/A * log(A/B) - (B - 1) * InvOnemBoA - 0.5 * (B + 1);
C *= 2.0 * InvOnemBoA; C *= 2.0 * InvOnemBoA;
k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value(); k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0;
partKsField_[cellI] = k / (1 - voidfraction); partKsField_[cellI] = k / (1 - voidfraction);
} }
} }

View File

@ -59,7 +59,32 @@ thermCondModel::thermCondModel
IOobject::NO_WRITE IOobject::NO_WRITE
) )
), ),
kf0_(transportProperties_.lookup("kf")), kf0Field_
(
IOobject
(
"kf0",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.)
),
CpField_
(
IOobject
(
"Cp",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.)
),
thermCondField_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("thCond"))), thermCondField_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("thCond"))),
wallQFactor_ wallQFactor_
( IOobject ( IOobject
@ -101,6 +126,38 @@ thermCondModel::thermCondModel
), ),
hasWallHeatLoss_(false) hasWallHeatLoss_(false)
{ {
// build constant fields for single phase case
if (!particleCloud_.multiphase())
{
kf0Field_ = volScalarField
(
IOobject
(
"kf0",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("kf"))
);
CpField_ = volScalarField
(
IOobject
(
"Cp",
particleCloud_.mesh().time().timeName(),
particleCloud_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
particleCloud_.mesh(),
dimensionedScalar(transportProperties_.lookup("Cp"))
);
}
if (wallQFactor_.headerOk()) if (wallQFactor_.headerOk())
{ {
Info << "Found field for scaling wall heat flux.\n" << endl; Info << "Found field for scaling wall heat flux.\n" << endl;
@ -179,6 +236,30 @@ void thermCondModel::calcBoundaryCorrections()
} }
} }
const volScalarField& thermCondModel::kf0Field() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("kf");
}
else
{
return kf0Field_;
}
}
const volScalarField& thermCondModel::CpField() const
{
if (particleCloud_.multiphase())
{
return particleCloud_.mesh().lookupObject<volScalarField>("Cp");
}
else
{
return CpField_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam } // End namespace Foam

View File

@ -63,7 +63,9 @@ protected:
IOdictionary transportProperties_; IOdictionary transportProperties_;
dimensionedScalar kf0_; mutable volScalarField kf0Field_;
mutable volScalarField CpField_;
volScalarField& thermCondField_; volScalarField& thermCondField_;
@ -132,6 +134,10 @@ public:
virtual void calcBoundaryCorrections(); virtual void calcBoundaryCorrections();
const volScalarField& kf0Field() const;
const volScalarField& CpField() const;
}; };

View File

@ -3,9 +3,11 @@ cfdemCloud = $(path)/cfdemCloud
derived = $(path)/derived derived = $(path)/derived
cfdTools = $(path)/cfdTools cfdTools = $(path)/cfdTools
energyModels = $(path)/subModels/energyModel energyModels = $(path)/subModels/energyModel
massTransferModels = $(path)/subModels/massTransferModel
forceModels = $(path)/subModels/forceModel forceModels = $(path)/subModels/forceModel
forceSubModels = $(path)/subModels/forceModel/forceSubModels forceSubModels = $(path)/subModels/forceModel/forceSubModels
thermCondModels = $(path)/subModels/thermCondModel thermCondModels = $(path)/subModels/thermCondModel
diffCoeffModels = $(path)/subModels/diffCoeffModel
chemistryModels = $(path)/subModels/chemistryModel chemistryModels = $(path)/subModels/chemistryModel
IOModels = $(path)/subModels/IOModel IOModels = $(path)/subModels/IOModel
voidFractionModels = $(path)/subModels/voidFractionModel voidFractionModels = $(path)/subModels/voidFractionModel
@ -42,14 +44,24 @@ $(energyModels)/energyModel/newEnergyModel.C
$(energyModels)/heatTransferGunn/heatTransferGunn.C $(energyModels)/heatTransferGunn/heatTransferGunn.C
$(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C $(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C
$(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C $(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C
$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C
$(energyModels)/reactionHeat/reactionHeat.C $(energyModels)/reactionHeat/reactionHeat.C
$(massTransferModels)/massTransferModel/massTransferModel.C
$(massTransferModels)/massTransferModel/newMassTransferModel.C
$(massTransferModels)/massTransferGunn/massTransferGunn.C
$(thermCondModels)/thermCondModel/thermCondModel.C $(thermCondModels)/thermCondModel/thermCondModel.C
$(thermCondModels)/thermCondModel/newThermCondModel.C $(thermCondModels)/thermCondModel/newThermCondModel.C
$(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C $(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
$(thermCondModels)/noTherm/noThermCond.C $(thermCondModels)/noTherm/noThermCond.C
$(diffCoeffModels)/diffCoeffModel/diffCoeffModel.C
$(diffCoeffModels)/diffCoeffModel/newDiffCoeffModel.C
$(diffCoeffModels)/SyamlalDiffCoeff/SyamlalDiffCoeff.C
$(diffCoeffModels)/noDiff/noDiffCoeff.C
$(forceModels)/forceModel/forceModel.C $(forceModels)/forceModel/forceModel.C
$(forceModels)/forceModel/newForceModel.C $(forceModels)/forceModel/newForceModel.C
$(forceModels)/noDrag/noDrag.C $(forceModels)/noDrag/noDrag.C
@ -65,6 +77,7 @@ $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
$(forceModels)/KochHillDrag/KochHillDrag.C $(forceModels)/KochHillDrag/KochHillDrag.C
$(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
$(forceModels)/virtualMassForce/virtualMassForce.C $(forceModels)/virtualMassForce/virtualMassForce.C
$(forceModels)/ParmarBassetForce/ParmarBassetForce.C
$(forceModels)/gradPForce/gradPForce.C $(forceModels)/gradPForce/gradPForce.C
$(forceModels)/viscForce/viscForce.C $(forceModels)/viscForce/viscForce.C
$(forceModels)/MeiLift/MeiLift.C $(forceModels)/MeiLift/MeiLift.C

View File

@ -0,0 +1,25 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
RASModel laminar;
turbulence off;
printCoeffs on;
// ************************************************************************* //

View File

@ -0,0 +1,304 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.4 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object couplingProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//===========================================================================//
// sub-models & settings
syncMode false;
//verbose;
useDDTvoidfraction;
modelType "A";
couplingInterval 5;
voidFractionModel divided;
locateModel engine;
meshMotionModel noMeshMotion;
regionModel allRegion;
IOModel basicIO;
probeModel off;
dataExchangeModel twoWayMPI;
averagingModel dense;
clockModel standardClock;
smoothingModel constDiffSmoothing;
getParticleDensities true;
getParticleAngVels true;
multiphase true;
forceModels
(
BeetstraDrag
MeiLift
virtualMassForce
ParmarBassetForce
viscForce
gradPForceSmooth
surfaceTensionForce
);
momCoupleModels
(
implicitCouple
);
energyModels
(
heatTransferGunn
wallHeatTransferYagi
);
massTransferModels
(
massTransferGunn
);
chemistryModels
(
off
);
thermCondModel SyamlalThermCond;
diffCoeffModel SyamlalDiffCoeff;
turbulenceModelType "turbulenceProperties";
//===========================================================================//
// sub-model properties
heatTransferGunnProps
{
interpolation true;
verbose false;
partRefTemp 1650;
calcPartTempField true;
partTempName "Temp";
partHeatFluxName "convectiveHeatFlux";
implicit true;
}
wallHeatTransferYagiProps
{
interpolation true;
verbose true;
wallTempName "wallTemp";
granVelFieldName "Us";
voidfractionMax 0.5;
implicit true;
}
massTransferGunnProps
{
interpolation false;
verbose false;
partMassFluxName "convectiveMassFlux";
implicit false;
coupleDEM true;
}
SyamlalThermCondProps
{
}
SyamlalDiffCoeffProps
{
}
BeetstraDragProps
{
velFieldName "U";
granVelFieldName "Us";
voidfractionFieldName "voidfraction";
interpolation true;
verbose false;
implForceDEM true;
}
dSauterProps
{
}
KochHillDragProps
{
velFieldName "U";
voidfractionFieldName "voidfraction";
implForceDEM true;
}
MeiLiftProps
{
velFieldName "U";
useShearInduced true;
useSpinInduced true;
combineShearSpin false;
interpolation true;
verbose false;
}
virtualMassForceProps
{
velFieldName "U";
voidfractionFieldName "voidfraction";
granVelFieldName "Us";
phiFieldName "phi";
useUs true;
useFelderhof true;
interpolation true;
smoothingModel constDiffAndTemporalSmoothing;
constDiffAndTemporalSmoothingProps
{
lowerLimit 1e-8;
upperLimit 1e8;
smoothingLength 0.05; // 2dp
smoothingStrength 0.0005; // timeScale = dt/gamma
correctBoundary true;
}
}
ParmarBassetForceProps
{
velFieldName "U";
granVelFieldName "Us";
useUs true;
interpolation true;
nIntegral 10;
discretisationOrder 1;
smoothingModel constDiffAndTemporalSmoothing;
constDiffAndTemporalSmoothingProps
{
lowerLimit 1e-8;
upperLimit 1e8;
smoothingLength 0.05; // 2dp
smoothingStrength 0.0005; // timeScale = dt/gamma
correctBoundary true;
}
}
viscForceProps
{
velocityFieldName "U";
}
gradPForceProps
{
pFieldName "p";
voidfractionFieldName "voidfraction";
velocityFieldName "U";
}
gradPForceSmoothProps
{
pFieldName "p_rgh";
voidfractionFieldName "voidfraction";
velocityFieldName "U";
smoothingModel "temporalSmoothing";
temporalSmoothingProps
{
lowerLimit 0.1;
upperLimit 1e10;
refField "p_rgh";
gamma 0.0005; // timeScale = dt/gamma
}
}
ArchimedesProps
{
gravityFieldName "g";
}
surfaceTensionForceProps
{
}
constDiffSmoothingProps
{
lowerLimit 0.1;
upperLimit 1e10;
smoothingLength 0.1; // 2dx
}
implicitCoupleProps
{
velFieldName "U";
granVelFieldName "Us";
voidfractionFieldName "voidfraction";
}
volWeightedAverageProps
{
scalarFieldNames
(
voidfraction
);
vectorFieldNames
(
);
upperThreshold 0.999;
lowerThreshold 0;
}
totalMomentumExchangeProps
{
implicitMomExFieldName "Ksl";
explicitMomExFieldName "none";
fluidVelFieldName "U";
granVelFieldName "Us";
}
particleCellVolumeProps
{
upperThreshold 0.999;
lowerThreshold 0.;
}
engineProps
{
treeSearch true;
}
dividedProps
{
alphaMin 0.01;
//porosity 0.5;
//interpolation true;
}
twoWayMPIProps
{
liggghtsPath "../DEM/in.liggghts_run";
}
trilinearProps
{
alphaMin 0.01;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,22 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class uniformDimensionedVectorField;
location "constant";
object g;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -2 0 0 0 0];
value ( 0 0 -9.81 );
// ************************************************************************* //

View File

@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.4 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object liggghtsCommands;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
liggghtsCommandModels
(
execute
execute
execute
execute
runLiggghts
writeLiggghts
);
// ************************************************************************* //
/*runLiggghtsProps
{
preNo false;
}*/
writeLiggghtsProps
{
writeLast off;
writeName "post/restart/liggghts.restartCFDEM";
overwrite on;
}
executeProps0
{
command
(
delete_atoms
porosity
raceways
noBlanks
dollar
curlyOpen
fDelRun
curlyClose
blanks
label
compress
no
);
labels
(
123457
);
runFirst off;
runLast off;
runEveryCouplingStep off;
runEveryWriteStep off;
startTime 0.0;
endTime 10000;
timeInterval 0.1;
}
executeProps1
{
command
(
write_restart
noBlanks
dotdot
slash
DEM
slash
post
slash
restart
slash
liggghts.restartSequence
);
runFirst off;
runLast on;
runEveryCouplingStep off;
runEveryWriteStep off;
}
executeProps2
{
command
(
delete_atoms
region
freeboard
compress
no
);
runFirst off;
runLast off;
runEveryCouplingStep on;
runEveryWriteStep off;
}
executeProps3
{
command
(
set
region
raceways
noBlanks
property
slash
atom
blanks
Temp
noBlanks
dollar
curlyOpen
Tinlet
curlyClose
);
runFirst off;
runLast off;
runEveryCouplingStep on;
runEveryWriteStep off;
}

View File

@ -0,0 +1,51 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//Properties from https://www.jstage.jst.go.jp/article/tetsutohagane/100/8/100_925/_pdf
phases
(
iron
{
transportModel Newtonian;
nu nu [ 0 2 -1 0 0 0 0 ] 1e-06;
rho rho [ 1 -3 0 0 0 0 0 ] 6700;
Cp Cp [ 0 2 -2 -1 0 0 0 ] 850;
kf kf [ 1 1 -3 -1 0 0 0 ] 330; // sped up! 16.5*20;
D D [ 0 2 -1 0 0 0 0 ] 0.17e-6; // sped up! 0.85e-8*20;
Cs Cs [ 1 -3 0 0 0 0 0 ] 335;
}
air
{
transportModel Newtonian;
nu nu [ 0 2 -1 0 0 0 0 ] 3e-04;
rho rho [ 1 -3 0 0 0 0 0 ] 0.25;
Cp Cp [ 0 2 -2 -1 0 0 0 ] 1200;
kf kf [ 1 1 -3 -1 0 0 0 ] 2.0; // sped up! 0.1*20;
D D [ 0 2 -1 0 0 0 0 ] 0.17e-6; // sped up! 0.85e-8*20; // dummy
Cs Cs [ 1 -3 0 0 0 0 0 ] 0.00025;
}
);
sigmas
(
(iron air) 1.65 // N/m
);
// ************************************************************************* //

View File

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "1";
object C;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -3 0 0 0 0 0];
internalField uniform 0; // see setFieldsDict
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 135;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object Ksl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -3 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "1";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 1623.15;
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 1823.15;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
atmosphere
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}
walls
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type zeroGradient;
}
inlet
{
type flowRateInletVelocity;
volumetricFlowRate 0.00037; // (1 kg coke/s) / (400 kg coke/tonne iron) / (6.7 tonne/m3) = 3.7e-4 m3/s
value uniform (0 0 0); // placeholder
}
}
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object Us;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object alpha.air;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type fixedValue;
value uniform 1;
}
walls
{
type alphaContactAngle;
thetaProperties 1 (( iron air ) 90 0 0 0);
value uniform 0;
}
outlet
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,45 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object alpha.iron;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type fixedValue;
value uniform 0;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 1;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object pSmoothField;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0]; //Pa
internalField uniform 0;
boundaryField
{
walls
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
atmosphere
{
type fixedValue;
value uniform 140000;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0]; //Pa
internalField uniform 0;
boundaryField
{
atmosphere
{
type fixedValue;
value uniform 140000;
}
walls
{
type fixedFluxPressure;
value uniform 0;
}
outlet
{
type fixedValue;
value uniform 0;
}
inlet
{
type fixedFluxPressure;
value uniform 0;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object sSmoothField;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object vSmoothField;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 0 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,45 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object voidfraction;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type fixedValue;
value uniform 1;
}
walls
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 1;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object wallTemp;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type zeroGradient;
}
walls
{
type fixedValue;
value uniform 1473.15;
}
outlet
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,129 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application pisoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 60.0;
deltaT 0.0005;
writeControl adjustableRunTime;
writeInterval 0.5;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep no;
maxCo 0.1;
maxAlphaCo 0.1;
libs
(
);
functions
{
#includeFunc flowRatePatch(name=inlet)
#includeFunc flowRatePatch(name=outlet)
#includeFunc patchAverage(alpha.iron,name=inlet)
#includeFunc patchAverage(alpha.iron,name=outlet)
liquidLevel
{
functionObjectLibs ("libutilityFunctionObjects.so");
type coded;
name integral;
executeControl writeTime;
executeInterval 1;
writeControl writeTime;
writeInterval 1;
codeExecute
#{
scalar hearthRadius = 0.5;
scalar hearthArea = M_PI*hearthRadius*hearthRadius;
const volScalarField& alphaIron = mesh().lookupObject<volScalarField>("alpha.iron");
scalar liquidVol = 0;
forAll (mesh().V(), cellI)
{
liquidVol += alphaIron[cellI]*mesh().V()[cellI];
}
reduce(liquidVol, sumOp<scalar>());
scalar liquidLevel = liquidVol/hearthArea;
Info << "liquid level: " << liquidLevel << " m" << endl;
if(Pstream::master()) {
// output file
if(!isDir(mesh().time().path()/".."/"postProcessing"))
mkDir(mesh().time().path()/".."/"postProcessing");
fileName outputFile(mesh().time().path()/".."/"postProcessing"/"liquidLevel.dat");
std::ofstream file;
// header
if(!isFile(outputFile)) {
file.open(outputFile,ios_base::out);
file << "#time(s) \t liquid_lev(m)" << std::endl;
file.close();
}
// output to file
file.open(outputFile,ios_base::out | ios_base::app);
file << mesh().time().value() << " \t";
file << liquidLevel << std::endl;
file.close();
}
#};
}
}
// ************************************************************************* //

Some files were not shown because too many files have changed in this diff Show More