diff --git a/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C index a1ea1537..713a6b36 100644 --- a/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C @@ -55,14 +55,21 @@ Foam::multiphaseMixture::calcNu() const { PtrDictionary::const_iterator iter = phases_.begin(); + // 1/nu + tmp tnuInv = iter()/iter().nu(); + volScalarField& nuInv = tnuInv.ref(); + + // nu tmp tnu = iter()*iter().nu(); volScalarField& nu = tnu.ref(); for (++iter; iter != phases_.end(); ++iter) { - nu += iter()*iter().nu(); + nuInv += iter()/iter().nu(); } + nu = 1/nuInv; + return tnu; } diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/Allwclean b/applications/solvers/cfdemSolverMultiphaseScalar/Allwclean new file mode 100755 index 00000000..1c8a5a68 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory +set -x + +wclean libso multiphaseMixture +wclean + +#------------------------------------------------------------------------------ diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/Allwmake b/applications/solvers/cfdemSolverMultiphaseScalar/Allwmake new file mode 100755 index 00000000..fbe71a59 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/Allwmake @@ -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 + +#------------------------------------------------------------------------------ diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H b/applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H new file mode 100644 index 00000000..b9392817 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H @@ -0,0 +1,26 @@ + // get mixture properties + Cp = mixture.Cp(); + kf = mixture.kf(); + keff=particleCloud.thermCondM().thermCond(); + + // get scalar source from DEM + particleCloud.energyContributions(Qsource); + particleCloud.energyCoefficients(QCoeff); + + fvScalarMatrix EEqn + ( + rho*Cp*(fvm::ddt(voidfraction,T) + + fvm::div(phi,T)) + - fvc::div(keff*voidfraction*fvc::grad(T)) + == + Qsource + fvm::Sp(QCoeff,T) + ); + + + EEqn.relax(); + fvOptions.constrain(EEqn); + EEqn.solve(); + + particleCloud.clockM().start(31,"postFlow"); + particleCloud.postFlow(); + particleCloud.clockM().stop("postFlow"); diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/Make/files b/applications/solvers/cfdemSolverMultiphaseScalar/Make/files new file mode 100644 index 00000000..b4567eab --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/Make/files @@ -0,0 +1,3 @@ +cfdemSolverMultiphaseScalar.C + +EXE = $(CFDEM_APP_DIR)/cfdemSolverMultiphaseScalar diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/Make/options b/applications/solvers/cfdemSolverMultiphaseScalar/Make/options new file mode 100644 index 00000000..5d1853b8 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/Make/options @@ -0,0 +1,30 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +EXE_INC = \ + -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 \ + +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) diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/UEqn.H b/applications/solvers/cfdemSolverMultiphaseScalar/UEqn.H new file mode 100644 index 00000000..86d12e50 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/UEqn.H @@ -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); +} diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/additionalChecks.H b/applications/solvers/cfdemSolverMultiphaseScalar/additionalChecks.H new file mode 100644 index 00000000..6f485462 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/additionalChecks.H @@ -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; +} diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/alphaCourantNo.H b/applications/solvers/cfdemSolverMultiphaseScalar/alphaCourantNo.H new file mode 100644 index 00000000..61dafb14 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/alphaCourantNo.H @@ -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; + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C b/applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C new file mode 100644 index 00000000..0174eeee --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C @@ -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 . + + 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 "energyModel.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #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()); + 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" + + // --- 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; + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + particleCloud.clockM().stop("Flow"); + particleCloud.clockM().stop("Global"); + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/correctPhi.H b/applications/solvers/cfdemSolverMultiphaseScalar/correctPhi.H new file mode 100644 index 00000000..9afcd58a --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/correctPhi.H @@ -0,0 +1,11 @@ +CorrectPhi +( + U, + phi, + p_rgh, + dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1), + geometricZeroField(), + pimple +); + +#include "continuityErrs.H" diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/createFields.H b/applications/solvers/cfdemSolverMultiphaseScalar/createFields.H new file mode 100644 index 00000000..9a7d6281 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/createFields.H @@ -0,0 +1,246 @@ +//=============================== +// 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 keff +( + IOobject + ( + "keff", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0) +); + + +//======================== + +volScalarField rhoEps ("rhoEps", rho * voidfraction); + +// Construct incompressible turbulence model +autoPtr 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()); + + diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files new file mode 100644 index 00000000..998255e4 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files @@ -0,0 +1,5 @@ +phase/phase.C +alphaContactAngle/alphaContactAngleFvPatchScalarField.C +multiphaseMixture.C + +LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoamScalar diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options new file mode 100644 index 00000000..f8ffa1cf --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options @@ -0,0 +1,13 @@ +EXE_INC = \ + -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 + +LIB_LIBS = \ + -linterfaceProperties \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C new file mode 100644 index 00000000..a0d433f4 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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& iF +) +: + zeroGradientFvPatchScalarField(p, iF) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper), + thetaProps_(gcpsf.thetaProps_) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + zeroGradientFvPatchScalarField(p, iF), + thetaProps_(dict.lookup("thetaProperties")) +{ + evaluate(); +} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const DimensionedField& 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 + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H new file mode 100644 index 00000000..09249c72 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H @@ -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 . + +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& + ); + + //- Construct from patch, internal field and dictionary + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given alphaContactAngleFvPatchScalarField + // onto a new patch + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new alphaContactAngleFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + 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 + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.C new file mode 120000 index 00000000..b5cf3fd0 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.C @@ -0,0 +1 @@ +../alphaContactAngle/alphaContactAngleFvPatchScalarField.C \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.H new file mode 120000 index 00000000..c7e4030f --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/alphaContactAngleFvPatchScalarField.H @@ -0,0 +1 @@ +../alphaContactAngle/alphaContactAngleFvPatchScalarField.H \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.C new file mode 120000 index 00000000..139112c3 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.C @@ -0,0 +1 @@ +../multiphaseMixture.C \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.H new file mode 120000 index 00000000..d780b93e --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.H @@ -0,0 +1 @@ +../multiphaseMixture.H \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.C new file mode 120000 index 00000000..aba05d8b --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.C @@ -0,0 +1 @@ +../phase/phase.C \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.H new file mode 120000 index 00000000..0010ad9b --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/phase.H @@ -0,0 +1 @@ +../phase/phase.H \ No newline at end of file diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C new file mode 100644 index 00000000..faa3a1ec --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C @@ -0,0 +1,823 @@ +/*---------------------------------------------------------------------------*\ +License + + This is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This code is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with this code. If not, see . + + Copyright (C) 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, phases_, iter) + { + alphas_ += level*iter(); + level += 1.0; + } +} + + +Foam::tmp +Foam::multiphaseMixture::calcNu() const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + // 1/nu + tmp tnuInv = iter()/iter().nu(); + volScalarField& nuInv = tnuInv.ref(); + + // nu + tmp 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::multiphaseMixture::calcStf() const +{ + tmp 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, phases_, iter1) + { + const phase& alpha1 = iter1(); + + PtrDictionary::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::multiphaseMixture::rho() const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + tmp trho = iter()*iter().rho(); + volScalarField& rho = trho.ref(); + + for (++iter; iter != phases_.end(); ++iter) + { + rho += iter()*iter().rho(); + } + + return trho; +} + + +Foam::tmp +Foam::multiphaseMixture::rho(const label patchi) const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + tmp 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::multiphaseMixture::mu() const +{ + Info << "In multiphasemixture mu()" << endl; + return rho()*nu(); +// PtrDictionary::const_iterator iter = phases_.begin(); + +// tmp 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::multiphaseMixture::mu(const label patchi) const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + tmp 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::multiphaseMixture::muf() const +{ + + return nuf()*fvc::interpolate(rho()); +// PtrDictionary::const_iterator iter = phases_.begin(); + +// tmp 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::multiphaseMixture::nu() const +{ + return nu_; +} + +Foam::tmp +Foam::multiphaseMixture::nu(const label patchi) const +{ + //return nu_.boundaryField()[patchi]; + PtrDictionary::const_iterator iter = phases_.begin(); + + tmp 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::multiphaseMixture::nuf() const +{ + //return muf()/fvc::interpolate(rho()); + PtrDictionary::const_iterator iter = phases_.begin(); + + tmp 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::multiphaseMixture::Cp() const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + // rho*Cp + tmp trhoCp = iter()*iter().Cp()*iter().rho(); + volScalarField& rhoCp = trhoCp.ref(); + + // Cp + tmp 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::multiphaseMixture::kf() const +{ + PtrDictionary::const_iterator iter = phases_.begin(); + + // rho*Cp/kf + tmp trhoCpkf = iter()*iter().rho()*iter().Cp()/iter().kf(); + volScalarField& rhoCpkf = trhoCpkf.ref(); + + // kf + tmp 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; +} + +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 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, phases_, iter) + { + iter().correct(); + } +} + + +Foam::tmp 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::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(gbf[patchi])) + { + const alphaContactAngleFvPatchScalarField& acap = + refCast(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::multiphaseMixture::K +( + const phase& alpha1, + const phase& alpha2 +) const +{ + tmp tnHatfv = nHatfv(alpha1, alpha2); + + correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef()); + + // Simple expression for curvature + return -fvc::div(tnHatfv & mesh_.Sf()); +} + + +Foam::tmp +Foam::multiphaseMixture::nearInterface() const +{ + tmp tnearInt + ( + new volScalarField + ( + IOobject + ( + "nearInterface", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("nearInterface", dimless, 0.0) + ) + ); + + forAllConstIter(PtrDictionary, 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 alphaPhiCorrs(phases_.size()); + int phasei = 0; + + forAllIter(PtrDictionary, 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, 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(), + 1, + 0, + 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, phases_, iter) + { + phase& alpha = iter(); + + surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei]; + alphaPhi += upwind(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 phaseData(lookup("phases")); + label phasei = 0; + + forAllIter(PtrDictionary, phases_, iter) + { + readOK &= iter().read(phaseData[phasei++].dict()); + } + + lookup("sigmas") >> sigmas_; + + return readOK; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H new file mode 100644 index 00000000..c7334d4a --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H @@ -0,0 +1,290 @@ +/*---------------------------------------------------------------------------*\ +License + + This is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This code is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with this code. If not, see . + + Copyright (C) 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 + { + public: + + class hash + : + public Hash + { + 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(alpha1Name, alpha2Name) + {} + + interfacePair(const phase& alpha1, const phase& alpha2) + : + Pair(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 phases_; + + const fvMesh& mesh_; + const volVectorField& U_; + const surfaceScalarField& phi_; + const volScalarField& voidfraction_; + surfaceScalarField rhoPhi_; + surfaceScalarField surfaceTensionForce_; + volScalarField alphas_; + + volScalarField nu_; + + typedef HashTable + 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 calcNu() const; + + void solveAlphas(const scalar cAlpha); + + tmp nHatfv + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + tmp nHatf + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + void correctContactAngle + ( + const phase& alpha1, + const phase& alpha2, + surfaceVectorField::Boundary& nHatb + ) const; + + tmp K(const phase& alpha1, const phase& alpha2) const; + tmp 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& 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 rho() const; + + //- Return the mixture density for patch + tmp rho(const label patchi) const; + + //- Return the dynamic laminar viscosity + tmp mu() const; + + //- Return the dynamic laminar viscosity for patch + tmp mu(const label patchi) const; + + //- Return the face-interpolated dynamic laminar viscosity + tmp muf() const; + + //- Return the kinematic laminar viscosity + tmp nu() const; + + //- Return the laminar viscosity for patch + tmp nu(const label patchi) const; + + //- Return the face-interpolated dynamic laminar viscosity + tmp nuf() const; + + //- Return the heat capacity + tmp Cp() const; + + //- Return the thermal conductivity + tmp kf() const; + + tmp surfaceTensionForce() const + { + return surfaceTensionForce_; + } + + //- Indicator of the proximity of the interface + // Field values are 1 near and 0 away for the interface. + tmp 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 + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C new file mode 100644 index 00000000..5b0d90ec --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#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_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::autoPtr Foam::phase::clone() const +{ + NotImplemented; + return autoPtr(NULL); +} + + +void Foam::phase::correct() +{ + nuModel_->correct(); +} + + +bool Foam::phase::read(const dictionary& phaseDict) +{ + phaseDict_ = phaseDict; + + phaseDict_.lookup("Cp") >> Cp_; + phaseDict_.lookup("kf") >> kf_; + + if (nuModel_->read(phaseDict_)) + { + phaseDict_.lookup("rho") >> rho_; + + return true; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H new file mode 100644 index 00000000..e61f63c9 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H @@ -0,0 +1,176 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +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 nuModel_; + dimensionedScalar rho_; + dimensionedScalar Cp_; + dimensionedScalar kf_; + +public: + + // Constructors + + //- Construct from components + phase + ( + const word& name, + const dictionary& phaseDict, + const volVectorField& U, + const surfaceScalarField& phi + ); + + //- Return clone + autoPtr 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 operator()(Istream& is) const + { + dictionaryEntry ent(dictionary::null, is); + return autoPtr(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 nu() const + { + return nuModel_->nu(); + } + + //- Return the laminar viscosity for patch + tmp 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_; + } + + //- 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 + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphaseScalar/pEqn.H b/applications/solvers/cfdemSolverMultiphaseScalar/pEqn.H new file mode 100644 index 00000000..e5a374f0 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphaseScalar/pEqn.H @@ -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; + } +} diff --git a/applications/solvers/cfdemSolverPiso/pEqn.H b/applications/solvers/cfdemSolverPiso/pEqn.H index c50eb9e3..8cce9907 100644 --- a/applications/solvers/cfdemSolverPiso/pEqn.H +++ b/applications/solvers/cfdemSolverPiso/pEqn.H @@ -21,7 +21,7 @@ phi = phiHbyA + rAUf*(fvc::interpolate(Ksl/rho) * phiS); MRF.makeRelative(phi); // adjustment of phi (only in cases w.o. p boundary conditions) not tested yet -adjustPhi(phi, U, p); +//adjustPhi(phi, U, p); // Update the pressure BCs to ensure flux consistency constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF); diff --git a/etc/bashrc b/etc/bashrc index 97057592..3d8c6581 100755 --- a/etc/bashrc +++ b/etc/bashrc @@ -218,7 +218,7 @@ alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh' #- 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 alias touchRec='find ./* -exec touch {} \;' diff --git a/etc/solver-list.txt b/etc/solver-list.txt index c370a4d4..54da1bb5 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -6,3 +6,4 @@ cfdemSolverIB/dir cfdemSolverPisoScalar/dir cfdemSolverRhoPimpleChem/dir cfdemSolverMultiphase/dir +cfdemSolverMultiphaseScalar/dir diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C index 00f0033b..0b96120e 100755 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C @@ -48,7 +48,8 @@ uniformFixedValueTubeFvPatchField::uniformFixedValueTubeFvPatchField densityFieldName_("rho"), tubeLength_(-1.), tubeDiameter_(-1.), - zeta_(-1) + zeta_(-1), + p0_(-1) { Info << "error, wrong constructor1!" << endl; } @@ -71,7 +72,8 @@ uniformFixedValueTubeFvPatchField::uniformFixedValueTubeFvPatchField densityFieldName_("rho"), tubeLength_(ptf.tubeLength_), tubeDiameter_(ptf.tubeDiameter_), - zeta_(ptf.zeta_) + zeta_(ptf.zeta_), + p0_(ptf.p0_) { const scalar t = this->db().time().timeOutputValue(); fvPatchField::operator==(uniformValue_->value(t)); @@ -89,13 +91,14 @@ uniformFixedValueTubeFvPatchField::uniformFixedValueTubeFvPatchField : fixedValueFvPatchField(p, iF), uniformValue_(Function1::New("uniformValue", dict)), - pName_("p"), //JOKER + pName_("p_rgh"), //JOKER phiName_("phi"), //JOKER velocityFieldName_("U"), densityFieldName_("rho"), tubeLength_(readScalar(dict.lookup("tubeLength"))), 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(); fvPatchField::operator==(uniformValue_->value(t)); @@ -110,13 +113,14 @@ uniformFixedValueTubeFvPatchField::uniformFixedValueTubeFvPatchField : fixedValueFvPatchField(ptf), uniformValue_(ptf.uniformValue_().clone().ptr()), - pName_("p"), //JOKER + pName_("p_rgh"), //JOKER phiName_("phi"), //JOKER velocityFieldName_("U"), densityFieldName_("rho"), tubeLength_(ptf.tubeLength_), tubeDiameter_(ptf.tubeDiameter_), - zeta_(ptf.zeta_) + zeta_(ptf.zeta_), + p0_(ptf.p0_) { const scalar t = this->db().time().timeOutputValue(); fvPatchField::operator==(uniformValue_->value(t)); @@ -132,13 +136,14 @@ uniformFixedValueTubeFvPatchField::uniformFixedValueTubeFvPatchField : fixedValueFvPatchField(ptf, iF), uniformValue_(ptf.uniformValue_().clone().ptr()), - pName_("p"), //JOKER + pName_("p_rgh"), //JOKER phiName_("phi"), //JOKER velocityFieldName_("U"), densityFieldName_("rho"), tubeLength_(ptf.tubeLength_), tubeDiameter_(ptf.tubeDiameter_), - zeta_(ptf.zeta_) + zeta_(ptf.zeta_), + p0_(ptf.p0_) { const scalar t = this->db().time().timeOutputValue(); fvPatchField::operator==(uniformValue_->value(t)); @@ -181,11 +186,11 @@ void uniformFixedValueTubeFvPatchField::updateCoeffs() //calc cell velocity from flux phip //scalar vel = phip/this->patch().magSf(); // some relaxation might be useful? - + // calc pressure drop // fvPatchField::operator==(0.5*zeta*density*mag(velocity)*mag(velocity)*uniformValue_->value(t)); //dP = zeta * l/d * rho u^2 / 2 - fvPatchField::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()-pressure)*uRelFact)*uniformValue_->value(t)); + fvPatchField::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()+p0_-pressure)*uRelFact)*uniformValue_->value(t)); //fvPatchField::operator==((100.*uniformValue_->value(t))); fixedValueFvPatchField::updateCoeffs(); } @@ -200,6 +205,7 @@ void uniformFixedValueTubeFvPatchField::write(Ostream& os) const os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl; os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl; os.writeKeyword("zeta") << zeta_ << token::END_STATEMENT << nl; + os.writeKeyword("p0") << p0_ << token::END_STATEMENT << nl; } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H index d86a92f2..ab4f859f 100755 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H @@ -74,6 +74,8 @@ class uniformFixedValueTubeFvPatchField scalar zeta_; + scalar p0_; + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 2d409055..aad43ee7 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -41,6 +41,8 @@ $(energyModels)/energyModel/energyModel.C $(energyModels)/energyModel/newEnergyModel.C $(energyModels)/heatTransferGunn/heatTransferGunn.C $(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C +$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C +$(energyModels)/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C $(energyModels)/reactionHeat/reactionHeat.C $(thermCondModels)/thermCondModel/thermCondModel.C @@ -65,6 +67,7 @@ $(forceModels)/KochHillDrag/KochHillDrag.C $(forceModels)/KochHillRWDrag/KochHillRWDrag.C $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C $(forceModels)/virtualMassForce/virtualMassForce.C +$(forceModels)/ParmarBassetForce/ParmarBassetForce.C $(forceModels)/gradPForce/gradPForce.C $(forceModels)/viscForce/viscForce.C $(forceModels)/MeiLift/MeiLift.C @@ -174,5 +177,6 @@ $(smoothingModels)/smoothingModel/newSmoothingModel.C $(smoothingModels)/noSmoothing/noSmoothing.C $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C $(smoothingModels)/temporalSmoothing/temporalSmoothing.C +$(smoothingModels)/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index bb205cc4..fbc4e304 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -87,6 +87,7 @@ cfdemCloud::cfdemCloud getParticleDensities_(couplingProperties_.lookupOrDefault("getParticleDensities",false)), getParticleEffVolFactors_(couplingProperties_.lookupOrDefault("getParticleEffVolFactors",false)), getParticleTypes_(couplingProperties_.lookupOrDefault("getParticleTypes",false)), + getParticleAngVels_(couplingProperties_.lookupOrDefault("getParticleAngVels",false)), maxDEMForce_(0.), modelType_(couplingProperties_.lookup("modelType")), positions_(NULL), @@ -103,6 +104,7 @@ cfdemCloud::cfdemCloud particleDensities_(NULL), particleEffVolFactors_(NULL), particleTypes_(NULL), + particleAngVels_(NULL), particleWeights_(NULL), particleVolumes_(NULL), particleV_(NULL), @@ -390,6 +392,7 @@ cfdemCloud::~cfdemCloud() if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1); if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1); if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1); + if(getParticleAngVels_) dataExchangeM().destroy(particleAngVels_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -405,6 +408,7 @@ void cfdemCloud::getDEMdata() if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_); if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_); if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_); + if(getParticleAngVels_) dataExchangeM().getData("omega","vector-atom",particleAngVels_); } void cfdemCloud::giveDEMdata() @@ -448,7 +452,12 @@ void cfdemCloud::setForces() resetArray(expForces_,numberOfParticles(),3); resetArray(DEMForces_,numberOfParticles(),3); resetArray(Cds_,numberOfParticles(),1); - for (int i=0;i("multiphase",false)), + kfFieldName_(propsDict_.lookupOrDefault("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + kfField_(sm.mesh().lookupObject (kfFieldName_)), + CpFieldName_(propsDict_.lookupOrDefault("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + CpField_(sm.mesh().lookupObject (CpFieldName_)) { allocateMyArrays(); @@ -341,7 +346,15 @@ void heatTransferGunn::calcEnergyContribution() ds = 2.*particleCloud_.radius(index); ds_scaled = ds/cg; muf = mufField[cellI]; + Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf; + + if(multiphase_) + { + kf0_ = kfField_[cellI]; + Cp_ = CpField_[cellI]; + } + Pr = max(SMALL, Cp_ * muf / kf0_); Nup = Nusselt(voidfraction, Rep, Pr); diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H index 3a03a480..71c5a9c9 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H @@ -116,6 +116,16 @@ protected: scalarList typeCG_; + bool multiphase_; + + word kfFieldName_; + + const volScalarField& kfField_; + + word CpFieldName_; + + const volScalarField& CpField_; + void allocateMyArrays() const; void partTempField(); diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.C b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.C new file mode 100644 index 00000000..b7c93ee4 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.C @@ -0,0 +1,329 @@ +/*---------------------------------------------------------------------------*\ +License + +This is free software: you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This code is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this code. If not, see . + +Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria + +\*---------------------------------------------------------------------------*/ + +#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")), + interpolation_(propsDict_.lookupOrDefault("interpolation",false)), + verbose_(propsDict_.lookupOrDefault("verbose",false)), + QWallFluidName_(propsDict_.lookupOrDefault("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) + ), + 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) + ), + 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("tempFieldName","T")), + tempField_(sm.mesh().lookupObject (tempFieldName_)), + voidfractionFieldName_(propsDict_.lookupOrDefault("voidfractionFieldName","voidfraction")), + voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + voidfractionMax_(readScalar(propsDict_.lookup("voidfractionMax"))), + maxSource_(1e30), + velFieldName_(propsDict_.lookupOrDefault("velFieldName","U")), + U_(sm.mesh().lookupObject (velFieldName_)), + UsFieldName_(propsDict_.lookup("granVelFieldName")), + Us_(sm.mesh().lookupObject (UsFieldName_)), + densityFieldName_(propsDict_.lookupOrDefault("densityFieldName","rho")), + rho_(sm.mesh().lookupObject (densityFieldName_)), + partRe_(NULL), + multiphase_(propsDict_.lookupOrDefault("multiphase",false)), + kfFieldName_(propsDict_.lookupOrDefault("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + kfField_(sm.mesh().lookupObject (kfFieldName_)), + CpFieldName_(propsDict_.lookupOrDefault("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + CpField_(sm.mesh().lookupObject (CpFieldName_)) + { + allocateMyArrays(); + + if (propsDict_.found("maxSource")) + { + maxSource_=readScalar(propsDict_.lookup ("maxSource")); + Info << "limiting wall source field to: " << maxSource_ << endl; + } + + if (verbose_) + { + ReField_.writeOpt() = IOobject::AUTO_WRITE; + PrField_.writeOpt() = IOobject::AUTO_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() + { + particleCloud_.dataExchangeM().destroy(partRe_,1); + } + + // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // + void wallHeatTransferYagi::allocateMyArrays() const + { + // get memory for 2d arrays + double initVal=0.0; + + if(verbose_) + { + particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1); + } + } + + // * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // + + void wallHeatTransferYagi::calcEnergyContribution() + { + allocateMyArrays(); + + // reset Scalar field + QWallFluid_.primitiveFieldRef() = 0.0; + + #ifdef compre + const volScalarField mufField = particleCloud_.turbulence().mu(); + #else + const volScalarField mufField = particleCloud_.turbulence().nu()*rho_; + #endif + + interpolationCellPoint voidfractionInterpolator_(voidfraction_); + interpolationCellPoint UInterpolator_(U_); + interpolationCellPoint TInterpolator_(tempField_); + + // calculate Rep + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + label cellI = particleCloud_.cellIDs()[index][0]; + if(cellI >= 0) + { + scalar voidfraction; + vector Ufluid; + + if(interpolation_) + { + vector position = particleCloud_.position(index); + voidfraction = voidfractionInterpolator_.interpolate(position,cellI); + Ufluid = UInterpolator_.interpolate(position,cellI); + } + else + { + voidfraction = voidfraction_[cellI]; + Ufluid = U_[cellI]; + } + + if (voidfraction < 0.01) + voidfraction = 0.01; + + vector Us = particleCloud_.velocity(index); + scalar magUr = mag(Ufluid - Us); + scalar ds = 2.*particleCloud_.radius(index); + scalar muf = mufField[cellI]; + + scalar Rep = ds * magUr * voidfraction * rho_[cellI]/ muf; + partRe_[index][0] = Rep; + } + } + + // calculate Rep field + particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().setScalarAverage + ( + ReField_, + partRe_, + particleCloud_.particleWeights(), + particleCloud_.averagingM().UsWeightField(), + NULL + ); + + // calculate Pr field + PrField_ = CpField_ * mufField / kfField_; + + 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 Urel + scalar magG = mag(U_[faceCelli]-Us_[faceCelli])*voidfraction_[faceCelli]*rho_[faceCelli]; + + // calculate H + scalar H; + if (voidfraction_[faceCelli]<=voidfractionMax_) + H = 0.2087 * (pow(ReField_[faceCelli]+SMALL,-0.20)) * CpField_[faceCelli] * magG / (pow(PrField_[faceCelli],2/3) + SMALL); + else + H = 0; + + // get delta T (wall-fluid) + scalar Twall = wallTemp_.boundaryField()[patchi][facei]; + scalar Tfluid = tempField_[faceCelli]; + scalar deltaT = Twall - Tfluid; + + // get area + scalar area = curPatch.magSf()[facei]; + + // calculate heat flux + heatFlux(faceCelli, H, area, Twall, Tfluid); + heatFluxCoeff(faceCelli, H, area); + + if(verbose_ && facei >=0 && facei <2) + { + Info << "####################" << endl; + Info << "cellID: " << faceCelli << endl; + Info << "G : " << magG << endl; + Info << "Re: " << ReField_[faceCelli] << endl; + Info << "Pr: " << PrField_[faceCelli] << endl; + Info << "Cp: " << CpField_[faceCelli] << endl; + Info << "kf: " << kfField_[faceCelli] << endl; + Info << "H : " << H << endl; + Info << "Twall: " << Twall << endl; + Info << "Tfluid: " << Tfluid << endl; + Info << "dT: " << deltaT << endl; + Info << "q: " << H*deltaT << endl; + Info << "area: " << area << endl; + Info << "Q:" << H*deltaT*area << endl; + } + } + } + else + { + FatalError << "wallHeatTransferYagi requires zeroGradient BC for temperature field" << endl; + } + } + } + + QWallFluid_.primitiveFieldRef() /= QWallFluid_.mesh().V(); + + // limit source term + 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::heatFlux(label faceCelli, scalar H, scalar area, scalar Twall, scalar Tfluid) + { + QWallFluid_[faceCelli] += H * area * (Twall - Tfluid); + } + + void wallHeatTransferYagi::heatFluxCoeff(label faceCelli, scalar H, scalar area) + { + //no heat transfer coefficient in explicit model + } + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.H b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.H new file mode 100644 index 00000000..d3ea5e13 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ +License + + This is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This code is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with this code. If not, see . + + Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria + + 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 interpolation_; + + bool verbose_; + + word QWallFluidName_; + + volScalarField QWallFluid_; + + word wallTempName_; + + volScalarField wallTemp_; + + 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 UsFieldName_; + + const volVectorField& Us_; + + word densityFieldName_; + + const volScalarField& rho_; + + mutable double **partRe_; + + bool multiphase_; + + word kfFieldName_; + + const volScalarField& kfField_; + + word CpFieldName_; + + const volScalarField& CpField_; + + void allocateMyArrays() const; + + virtual void heatFlux(label, scalar, scalar, scalar, scalar); + + virtual void heatFluxCoeff(label, 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C new file mode 100644 index 00000000..9499fcba --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ +License + + This is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This code is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with this code. If not, see . + + Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "wallHeatTransferYagiImplicit.H" +#include "addToRunTimeSelectionTable.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(wallHeatTransferYagiImplicit, 0); + +addToRunTimeSelectionTable(energyModel, wallHeatTransferYagiImplicit, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +wallHeatTransferYagiImplicit::wallHeatTransferYagiImplicit +( + const dictionary& dict, + cfdemCloudEnergy& sm +) +: + wallHeatTransferYagi(dict,sm), + QWallFluidCoeffName_(propsDict_.lookupOrDefault("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) + ), + partHeatFluxCoeff_(NULL) +{ + allocateMyArrays(); + + // no limiting necessary for implicit heat transfer + maxSource_ = 1e30; +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +wallHeatTransferYagiImplicit::~wallHeatTransferYagiImplicit() +{ + particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1); +} + +// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // +void wallHeatTransferYagiImplicit::allocateMyArrays() const +{ +// wallHeatTransferYagi::allocateMyArrays(); + double initVal=0.0; + particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1); +} +// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // + +void wallHeatTransferYagiImplicit::calcEnergyContribution() +{ + allocateMyArrays(); + + QWallFluidCoeff_.primitiveFieldRef() = 0.0; + + wallHeatTransferYagi::calcEnergyContribution(); + + QWallFluidCoeff_.primitiveFieldRef() /= QWallFluidCoeff_.mesh().V(); + +// QWallFluidCoeff_.correctBoundaryConditions(); + +} + +void wallHeatTransferYagiImplicit::addEnergyCoefficient(volScalarField& Qsource) const +{ + Qsource += QWallFluidCoeff_; +} + +void wallHeatTransferYagiImplicit::heatFlux(label faceCelli, scalar H, scalar area, scalar Twall, scalar Tfluid) +{ + QWallFluid_[faceCelli] += H * area * Twall; +} + +void wallHeatTransferYagiImplicit::heatFluxCoeff(label faceCelli, scalar H, scalar area) +{ + QWallFluidCoeff_[faceCelli] -= H * area; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.H b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.H new file mode 100644 index 00000000..08c8cd31 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.H @@ -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 . + + 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 wallHeatTransferYagiImplicit_H +#define wallHeatTransferYagiImplicit_H + +#include "fvCFD.H" +#include "cfdemCloudEnergy.H" +#include "wallHeatTransferYagi.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class wallHeatTransferYagiImplicit Declaration +\*---------------------------------------------------------------------------*/ + +class wallHeatTransferYagiImplicit +: + public wallHeatTransferYagi +{ +private: + + word QWallFluidCoeffName_; + + volScalarField QWallFluidCoeff_; + + mutable double **partHeatFluxCoeff_; // Lagrangian array + + void allocateMyArrays() const; + + virtual void heatFlux(label, scalar, scalar, scalar, scalar); + + virtual void heatFluxCoeff(label, scalar, scalar); + +public: + + //- Runtime type information + TypeName("wallHeatTransferYagiImplicit"); + + // Constructors + + //- Construct from components + wallHeatTransferYagiImplicit + ( + const dictionary& dict, + cfdemCloudEnergy& sm + ); + + + // Destructor + + virtual ~wallHeatTransferYagiImplicit(); + + + // Member Functions + + void addEnergyCoefficient(volScalarField&) const; + + void calcEnergyContribution(); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 18381aa3..01aea1e0 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -70,14 +70,21 @@ BeetstraDrag::BeetstraDrag k_(0.05), useGC_(false), usePC_(false) + polydisperse_(propsDict_.lookupOrDefault("polydisperse",false)) { + + if(polydisperse_) + Info << "Drag model: using polydisperse correction factor \n"; + + //Append the field names to be probed particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().scalarFields_.append("Rep"); - particleCloud_.probeM().scalarFields_.append("betaP"); particleCloud_.probeM().scalarFields_.append("voidfraction"); + particleCloud_.probeM().scalarFields_.append("Fdrag"); + particleCloud_.probeM().scalarFields_.append("y_polydisperse"); particleCloud_.probeM().writeHeader(); // init force sub model @@ -159,7 +166,7 @@ void BeetstraDrag::setForce() const scalar voidfraction(1); vector Ufluid(0,0,0); vector drag(0,0,0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -177,12 +184,14 @@ void BeetstraDrag::setForce() const scalar cg = typeCG_[0]; label partType = 1; + scalar yi(0); vector dragExplicit(0,0,0); scalar dragCoefficient(0); interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); + interpolationCellPoint dSauterInterpolator_(dSauterField_); #include "setupProbeModel.H" @@ -205,6 +214,7 @@ void BeetstraDrag::setForce() const Ufluid = UInterpolator_.interpolate(position,cellI); //Ensure interpolated void fraction to be meaningful // Info << " --> voidfraction: " << voidfraction << endl; + if (voidfraction > 1.00) voidfraction = 1.0; if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_; } @@ -223,6 +233,7 @@ void BeetstraDrag::setForce() const scaleDia3 = cg*cg*cg; } +<<<<<<< HEAD Us = particleCloud_.velocity(index); Ur = Ufluid-Us; magUr = mag(Ur); @@ -243,6 +254,13 @@ void BeetstraDrag::setForce() const *3*M_PI*nuf*rho*voidfraction *effDiameter(ds_scaled, cellI, index) *scaleDia3*scaleDrag_; + if (polydisperse_) + { + // Beetstra et. al (2007), Eq. (21) + yi = ds_scaled/dSauter; + dragCoefficient *= (voidfraction*yi + 0.064*voidfraction*yi*yi*yi); + } + // calculate filtering corrections if (useGC_) @@ -274,6 +292,7 @@ void BeetstraDrag::setForce() const { Pout << "cellI = " << cellI << endl; Pout << "index = " << index << endl; + Pout << "Ufluid = " << Ufluid << endl; Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; Pout << "ds = " << ds << endl; @@ -295,6 +314,8 @@ void BeetstraDrag::setForce() const vValues.append(Ur); sValues.append(Rep); sValues.append(voidfraction); + sValues.append(Fdrag); + sValues.append(yi); particleCloud_.probeM().writeProbe(index, sValues, vValues); } } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index cfeeea39..920c556d 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H @@ -102,7 +102,7 @@ protected: double h(double) const; - + mutable bool polydisperse_; public: diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C index c753e8c1..b5368c17 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C @@ -81,9 +81,29 @@ MeiLift::MeiLift propsDict_(dict.subDict(typeName + "Props")), velFieldName_(propsDict_.lookup("velFieldName")), U_(sm.mesh().lookupObject (velFieldName_)), - useSecondOrderTerms_(false) + useShearInduced_(propsDict_.lookupOrDefault("useShearInduced",true)), + useSpinInduced_(propsDict_.lookupOrDefault("useSpinInduced",false)), + combineShearSpin_(propsDict_.lookupOrDefault("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("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 setForceSubModels(propsDict_); @@ -100,10 +120,15 @@ MeiLift::MeiLift particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug - particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug - particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug - particleCloud_.probeM().scalarFields_.append("Rew"); //other are debug - particleCloud_.probeM().scalarFields_.append("J_star"); //other are debug + particleCloud_.probeM().vectorFields_.append("vorticity"); + particleCloud_.probeM().vectorFields_.append("Ang_velocity"); + particleCloud_.probeM().scalarFields_.append("Rep"); + 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(); } @@ -122,27 +147,44 @@ void MeiLift::setForce() const const volScalarField& nufField = forceSubM(0).nuField(); const volScalarField& rhoField = forceSubM(0).rhoField(); + // vectors vector position(0,0,0); vector lift(0,0,0); vector Us(0,0,0); vector Ur(0,0,0); + vector Omega(0,0,0); + vector vorticity(0,0,0); + + // properties scalar magUr(0); scalar magVorticity(0); + scalar magOmega(0); scalar ds(0); scalar nuf(0); scalar rho(0); scalar voidfraction(1); + + // dimensionless groups scalar Rep(0); scalar Rew(0); - scalar Cl(0); - scalar Cl_star(0); - scalar J_star(0); - scalar Omega_eq(0); - scalar alphaStar(0); - scalar epsilonSqr(0.0); - scalar epsilon(0); + + // shear induced 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_); interpolationCellPoint UInterpolator_(U_); @@ -152,14 +194,14 @@ void MeiLift::setForce() const for (int index = 0; index < particleCloud_.numberOfParticles(); ++index) { - //if(mask[index][0]) - //{ lift = vector::zero; label cellI = particleCloud_.cellIDs()[index][0]; if (cellI > -1) // particle Found { - Us = particleCloud_.velocity(index); + // properties + Us = particleCloud_.velocity(index); + Omega = particleCloud_.particleAngVel(index); if (forceSubM(0).interpolation()) { @@ -173,23 +215,23 @@ void MeiLift::setForce() const vorticity = vorticityField[cellI]; } - magUr = mag(Ur); - magVorticity = mag(vorticity); + ds = 2. * particleCloud_.radius(index); + nuf = nufField[cellI]; + rho = rhoField[cellI]; - if (magUr > 0 && magVorticity > 0) + magUr = mag(Ur); + Rep = ds*magUr/nuf; + + // shear-induced lift + if (useShearInduced_) { - ds = 2. * particleCloud_.radius(index); - nuf = nufField[cellI]; - rho = rhoField[cellI]; + magVorticity = mag(vorticity); + Rew = magVorticity*ds*ds/nuf; + omega_star = magVorticity * ds / magUr; + alphaStar = 0.5 * omega_star; + epsilonSqr = omega_star / Rep; + epsilon = sqrt(epsilonSqr); - // calc dimensionless properties - Rep = ds*magUr/nuf; - Rew = magVorticity*ds*ds/nuf; - - omega_star = magVorticity * ds / magUr; - alphaStar = 0.5 * omega_star; - epsilonSqr = omega_star / Rep; - epsilon = sqrt(epsilonSqr); //Basic model for the correction to the Saffman lift //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009) @@ -212,35 +254,76 @@ void MeiLift::setForce() const * (1.0 + tanh(2.5 * (log10(epsilon) + 0.191))) * (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) - if (useSecondOrderTerms_) + //Loth and Dorgan (2009), Eq (31), Eq (32) + Clshear = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model + } + + 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) { - 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; - } - - //Loth and Dorgan (2009), Eq (27) - lift = 0.125 * constant::mathematical::pi - * rho - * Cl - * magUr * Ur ^ vorticity / magVorticity - * ds * ds; - - if (modelType_ == "B") - { - voidfraction = particleCloud_.voidfraction(index); - lift /= voidfraction; + //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 + * rho + * Clcombined + * magUr * magUr + * (Ur ^ vorticity) / mag(Ur ^ vorticity) // force direction + * 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") + { + voidfraction = particleCloud_.voidfraction(index); + lift /= voidfraction; + } + //********************************** //SAMPLING AND VERBOSE OUTOUT @@ -250,6 +333,7 @@ void MeiLift::setForce() const Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; Pout << "vorticity = " << vorticity << endl; + Pout << "Omega = " << Omega << endl; Pout << "ds = " << ds << endl; Pout << "rho = " << rho << endl; Pout << "nuf = " << nuf << endl; @@ -258,6 +342,10 @@ void MeiLift::setForce() const Pout << "alphaStar = " << alphaStar << endl; Pout << "epsilon = " << epsilon << 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; } @@ -268,9 +356,14 @@ void MeiLift::setForce() const vValues.append(lift); //first entry must the be the force vValues.append(Ur); vValues.append(vorticity); + vValues.append(Omega); sValues.append(Rep); sValues.append(Rew); sValues.append(J_star); + sValues.append(Clshear); + sValues.append(Clspin_star); + sValues.append(Omega_eq); + sValues.append(Clcombined); particleCloud_.probeM().writeProbe(index, sValues, vValues); } // END OF SAMPLING AND VERBOSE OUTOUT @@ -279,7 +372,6 @@ void MeiLift::setForce() const } // write particle based data to global array forceSubM(0).partToArray(index,lift,vector::zero); - //} } } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H index 81ecfaee..4b43a0d0 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H @@ -83,7 +83,11 @@ private: const volVectorField& U_; - bool useSecondOrderTerms_; + bool useShearInduced_; + + bool useSpinInduced_; + + bool combineShearSpin_; public: diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.C new file mode 100644 index 00000000..9d6c310e --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.C @@ -0,0 +1,626 @@ +/*---------------------------------------------------------------------------*\ + 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 (velFieldName_)), + UsFieldName_(propsDict_.lookup("granVelFieldName")), + Us_(sm.mesh().lookupObject (UsFieldName_)), + nInt_(readScalar(propsDict_.lookup("nIntegral"))), + discOrder_(readScalar(propsDict_.lookup("discretisationOrder"))), + nHist_(nInt_+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(2*discOrder_+1,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(0,0,0)) + ), + 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 + ) + ) +{ + + // 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!" << endl; + + //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 UrelInterpolator_(Urel_); + interpolationCellPoint ddtUrelInterpolator_(ddtUrel_); + + // + interpolationCellPoint UrelNoSmoothInterpolator_(UrelNoSmooth_); + interpolationCellPoint 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<.25 || r>2.) + { + 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 (t00; 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(vector F1, 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=2*discOrder_; 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 ddtUrelHist_; + + mutable List rHist_; + + mutable List> FHist_; + + mutable double** gH0_; + + mutable double** tRef_; + + mutable double** mRef_; + + mutable double** lRef_; + + mutable volVectorField Urel_; + + mutable volVectorField ddtUrel_; + + autoPtr smoothingModel_; + + +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_; + } + + scalar calculateK0(scalar r, scalar dt) const; + + scalar trapWeight(int i, int n) const; + + void update_ddtUrelHist(vector ddtUrel, int index) const; + + void update_rHist(scalar r, int index) const; + + void update_FHist(vector F1, 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; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C index 122b50e7..926142b8 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C @@ -67,11 +67,37 @@ virtualMassForce::virtualMassForce propsDict_(dict.subDict(typeName + "Props")), velFieldName_(propsDict_.lookup("velFieldName")), U_(sm.mesh().lookupObject (velFieldName_)), + voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), + voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + UsFieldName_(propsDict_.lookup("granVelFieldName")), + Us_(sm.mesh().lookupObject (UsFieldName_)), phiFieldName_(propsDict_.lookup("phiFieldName")), phi_(sm.mesh().lookupObject (phiFieldName_)), UrelOld_(NULL), splitUrelCalculation_(propsDict_.lookupOrDefault("splitUrelCalculation",false)), - Cadd_(0.5) + useUs_(propsDict_.lookupOrDefault("useUs",false)), + useFelderhof_(propsDict_.lookupOrDefault("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 + ) + ) { if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0) @@ -98,17 +124,39 @@ virtualMassForce::virtualMassForce Cadd_ = readScalar(propsDict_.lookup("Cadd")); 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("getParticleDensities",false)) + FatalError << "Virtual mass model: useFelderhof=true requires getParticleDensities=true in couplingProperties" << abort(FatalError); + } particleCloud_.checkCG(true); //Append the field names to be probed particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force - particleCloud_.probeM().vectorFields_.append("Urel"); - particleCloud_.probeM().vectorFields_.append("UrelOld"); - particleCloud_.probeM().vectorFields_.append("ddtUrel"); + particleCloud_.probeM().vectorFields_.append("acceleration"); + // + particleCloud_.probeM().vectorFields_.append("accelerationNoSmooth"); + // + particleCloud_.probeM().scalarFields_.append("Cadd"); particleCloud_.probeM().scalarFields_.append("Vs"); particleCloud_.probeM().scalarFields_.append("rho"); + particleCloud_.probeM().scalarFields_.append("voidfraction"); + particleCloud_.probeM().scalarFields_.append("specificGravity"); particleCloud_.probeM().writeHeader(); } @@ -125,22 +173,28 @@ virtualMassForce::~virtualMassForce() void virtualMassForce::setForce() const { - reAllocArrays(); + if (!useUs_ || !splitUrelCalculation_ ) + reAllocArrays(); scalar dt = U_.mesh().time().deltaT().value(); - vector position(0,0,0); - 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()); + //acceleration field 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 - interpolationCellPoint UInterpolator_( U_); - interpolationCellPoint DDtUInterpolator_(DDtU_); + // + volVectorField DDtUrelNoSmooth_ = DDtUrel_; + interpolationCellPoint DDtUrelNoSmoothInterpolator_(DDtUrelNoSmooth_); + // + + // smoothen + smoothingM().smoothen(DDtUrel_); + + interpolationCellPoint UInterpolator_(U_); + interpolationCellPoint DDtUrelInterpolator_(DDtUrel_); + interpolationCellPoint voidfractionInterpolator_(voidfraction_); #include "setupProbeModel.H" @@ -149,62 +203,94 @@ void virtualMassForce::setForce() const for(int index = 0;index < particleCloud_.numberOfParticles(); index++) { 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]; if (cellI > -1) // particle Found { - + // particle position if(forceSubM(0).interpolation()) + position = particleCloud_.position(index); + + //********* acceleration value *********// + if(splitUrelCalculation_ || useUs_ ) { - position = particleCloud_.position(index); - Ufluid = UInterpolator_.interpolate(position,cellI); + // DDtUrel from acceleration field + if(forceSubM(0).interpolation()) + DDtUrel = DDtUrelInterpolator_.interpolate(position,cellI); + else + DDtUrel = DDtUrel_[cellI]; } else { - Ufluid = U_[cellI]; + // DDtUrel from UrelOld + vector Ufluid(0,0,0); + + // 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 + if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step + haveUrelOld_ = false; + else + haveUrelOld_ = true; + + vector UrelOld(0,0,0); + + for(int j=0;j<3;j++) + { + UrelOld[j] = UrelOld_[index][j]; + UrelOld_[index][j] = Urel[j]; + } + if(haveUrelOld_ ) //only compute force if we have old (relative) velocity + DDtUrel = (Urel-UrelOld)/dt; } - - if(splitUrelCalculation_) //if split, just use total derivative of fluid velocity - if(forceSubM(0).interpolation()) - { - DDtU = DDtUInterpolator_.interpolate(position,cellI); - } - else - { - DDtU = DDtU_[cellI]; - } - else - { - vector Us = particleCloud_.velocity(index); - Ur = Ufluid - Us; - } - - - //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 - haveUrelOld_ = false; - else - haveUrelOld_ = true; - - vector UrelOld(0.,0.,0.); - vector ddtUrel(0.,0.,0.); - for(int j=0;j<3;j++) - { - UrelOld[j] = UrelOld_[index][j]; - UrelOld_[index][j] = Ur[j]; - } - if(haveUrelOld_ ) //only compute force if we have old (relative) velocity - ddtUrel = (Ur-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 Vs = ds*ds*ds*M_PI/6; - scalar rho = forceSubM(0).rhoField()[cellI]; + + scalar Cadd; + + if (useFelderhof_) + { + scalar epsilons(0); + scalar logsg(0); - virtualMassForce = Cadd_ * rho * Vs * ddtUrel; + 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) { @@ -212,21 +298,33 @@ void virtualMassForce::setForce() const Pout << "position = " << particleCloud_.position(index) << endl; } - //Set value fields and write the probe + // Set value fields and write the probe if(probeIt_) { + // + vector DDtUrelNoSmooth(0,0,0); + if(forceSubM(0).interpolation()) + DDtUrelNoSmooth = DDtUrelInterpolator_.interpolate(position,cellI); + else + DDtUrelNoSmooth = DDtUrelNoSmooth_[cellI]; + // #include "setupProbeModelfields.H" vValues.append(virtualMassForce); //first entry must the be the force - vValues.append(Ur); - vValues.append(UrelOld); - vValues.append(ddtUrel); + vValues.append(DDtUrel); + // + vValues.append(DDtUrelNoSmooth); + // + sValues.append(Cadd); sValues.append(Vs); sValues.append(rho); + sValues.append(voidfraction); + sValues.append(sg); particleCloud_.probeM().writeProbe(index, sValues, vValues); } } else //particle not on this CPU - UrelOld_[index][0]=NOTONCPU; + if (!useUs_ || !splitUrelCalculation_ ) + UrelOld_[index][0] = NOTONCPU; // write particle based data to global array forceSubM(0).partToArray(index,virtualMassForce,vector::zero); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H index 28597d0b..a0dee63f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H @@ -67,6 +67,14 @@ private: const volVectorField& U_; + word voidfractionFieldName_; + + const volScalarField& voidfraction_; + + word UsFieldName_; + + const volVectorField& Us_; + word phiFieldName_; const surfaceScalarField& phi_; @@ -76,7 +84,15 @@ private: const bool splitUrelCalculation_; //indicator to split calculation of Urel between CFDEM and LIGGGHTS //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_; public: @@ -102,6 +118,11 @@ public: void setForce() const; void reAllocArrays() const; + + inline const smoothingModel& smoothingM() const + { + return smoothingModel_; + } }; diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C index 9c6236d8..86e15269 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C @@ -147,7 +147,7 @@ void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict) // if this makes troubles try floor((startTime_+SMALL)/.. as above firstCouplingStep_ = floor((startTime_+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 { @@ -279,6 +279,11 @@ void liggghtsCommandModel::parseCommandList(wordList& commandList,labelList& lab else if (add=="dot") add = "."; else if (add=="dotdot") 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 { add = ""; diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C new file mode 100644 index 00000000..a12e566f --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C @@ -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("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 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.H b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.H new file mode 100644 index 00000000..16d962de --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.H @@ -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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C index ffa82e4e..701476f6 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C @@ -67,7 +67,12 @@ SyamlalThermCond::SyamlalThermCond sm.mesh(), dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0) ), - hasWallQFactor_(false) + hasWallQFactor_(false), + multiphase_(propsDict_.lookupOrDefault("multiphase",false)), + kfFieldName_(propsDict_.lookupOrDefault("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + kfField_(sm.mesh().lookupObject (kfFieldName_)), + CpFieldName_(propsDict_.lookupOrDefault("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase + CpField_(sm.mesh().lookupObject (CpFieldName_)) { if (wallQFactor_.headerOk()) { @@ -106,12 +111,17 @@ tmp SyamlalThermCond::thermCond() const volScalarField& svf = tvf.ref(); - forAll(svf,cellI) { - if (1-voidfraction_[cellI] < SMALL) svf[cellI] = kf0_.value(); + scalar kf; + if (multiphase_) + kf = kfField_[cellI]; + else + kf = kf0_.value(); + + if (1-voidfraction_[cellI] < SMALL) svf[cellI] = kf; else if (voidfraction_[cellI] < SMALL) svf[cellI] = 0.0; - else svf[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf0_.value(); + else svf[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf; } // if a wallQFactor field is present, use it to scale heat transport through a patch @@ -127,7 +137,10 @@ tmp SyamlalThermCond::thermCond() const tmp SyamlalThermCond::thermDiff() const { - return thermCond()/(rho_*Cp_); + if (multiphase_) + return thermCond()/(rho_*CpField_); + else + return thermCond()/(rho_*Cp_); } diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H b/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H index 1d17d611..a20c4fee 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H @@ -68,6 +68,16 @@ private: bool hasWallQFactor_; + bool multiphase_; + + word kfFieldName_; + + const volScalarField& kfField_; + + word CpFieldName_; + + const volScalarField& CpField_; + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files index 554c42b1..5faa22fe 100644 --- a/src/lagrangian/cfdemParticleComp/Make/files +++ b/src/lagrangian/cfdemParticleComp/Make/files @@ -40,6 +40,8 @@ $(energyModels)/energyModel/energyModel.C $(energyModels)/energyModel/newEnergyModel.C $(energyModels)/heatTransferGunn/heatTransferGunn.C $(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C +$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C +$(energyModels)/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C $(energyModels)/reactionHeat/reactionHeat.C $(thermCondModels)/thermCondModel/thermCondModel.C @@ -63,6 +65,7 @@ $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C $(forceModels)/KochHillDrag/KochHillDrag.C $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C $(forceModels)/virtualMassForce/virtualMassForce.C +$(forceModels)/ParmarBassetForce/ParmarBassetForce.C $(forceModels)/gradPForce/gradPForce.C $(forceModels)/viscForce/viscForce.C $(forceModels)/MeiLift/MeiLift.C diff --git a/tutorials/cfdemSolverMultiphase/tankDrainage/CFD/system/fvSchemes b/tutorials/cfdemSolverMultiphase/tankDrainage/CFD/system/fvSchemes index d0a271ce..b8a0dc65 100644 --- a/tutorials/cfdemSolverMultiphase/tankDrainage/CFD/system/fvSchemes +++ b/tutorials/cfdemSolverMultiphase/tankDrainage/CFD/system/fvSchemes @@ -38,7 +38,7 @@ divSchemes div(((dev(grad(U).T())*rho)*dev(grad(U).T()))) Gauss linear; div((((((alpha.water*rho)*nu.water)|(alpha.water*rho))*rho)*dev(grad(U).T()))) Gauss linear; div(((nu*rho)*dev(grad(U).T()))) Gauss linear; - default Gauss linear; + default Gauss linear; } laplacianSchemes