diff --git a/.gitignore b/.gitignore index 1d1ae5ca..b4b35554 100644 --- a/.gitignore +++ b/.gitignore @@ -7,5 +7,6 @@ log.* *~ **/linux*Gcc*/ +**/.vscode lnInclude diff --git a/applications/solvers/cfdemSolverMultiphase/Allwclean b/applications/solvers/cfdemSolverMultiphase/Allwclean new file mode 100755 index 00000000..1c8a5a68 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/Allwmake b/applications/solvers/cfdemSolverMultiphase/Allwmake new file mode 100755 index 00000000..fbe71a59 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/Make/files b/applications/solvers/cfdemSolverMultiphase/Make/files new file mode 100644 index 00000000..2c615be1 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/Make/files @@ -0,0 +1,3 @@ +cfdemSolverMultiphase.C + +EXE = $(CFDEM_APP_DIR)/cfdemSolverMultiphase diff --git a/applications/solvers/cfdemSolverMultiphase/Make/options b/applications/solvers/cfdemSolverMultiphase/Make/options new file mode 100644 index 00000000..0cfccaff --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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)\ + -lcfdemMultiphaseInterFoam \ + -linterfaceProperties \ + -lincompressibleTransportModels \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling \ + -l$(CFDEM_LIB_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/cfdemSolverMultiphase/UEqn.H b/applications/solvers/cfdemSolverMultiphase/UEqn.H new file mode 100644 index 00000000..ace635f5 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/additionalChecks.H b/applications/solvers/cfdemSolverMultiphase/additionalChecks.H new file mode 100644 index 00000000..6f485462 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/alphaCourantNo.H b/applications/solvers/cfdemSolverMultiphase/alphaCourantNo.H new file mode 100644 index 00000000..61dafb14 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/cfdemSolverMultiphase.C b/applications/solvers/cfdemSolverMultiphase/cfdemSolverMultiphase.C new file mode 100644 index 00000000..42e514da --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/cfdemSolverMultiphase.C @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ +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 + cfdemSolverMultiphase + +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 "cfdemCloud.H" +#include "implicitCouple.H" +#include "clockModel.H" +#include "smoothingModel.H" +#include "forceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "correctPhi.H" + #include "CourantNo.H" + + turbulence->validate(); + + // create cfdemCloud + cfdemCloud 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; + + // --- 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/cfdemSolverMultiphase/correctPhi.H b/applications/solvers/cfdemSolverMultiphase/correctPhi.H new file mode 100644 index 00000000..9afcd58a --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/createFields.H b/applications/solvers/cfdemSolverMultiphase/createFields.H new file mode 100644 index 00000000..3b4194c5 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/createFields.H @@ -0,0 +1,156 @@ +//=============================== +// 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(); + +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/cfdemSolverMultiphase/multiphaseMixture/Make/files b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/files new file mode 100644 index 00000000..572171ec --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/files @@ -0,0 +1,5 @@ +phase/phase.C +alphaContactAngle/alphaContactAngleFvPatchScalarField.C +multiphaseMixture.C + +LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoam diff --git a/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/options b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/options new file mode 100644 index 00000000..f8ffa1cf --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C new file mode 100644 index 00000000..a0d433f4 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H new file mode 100644 index 00000000..09249c72 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C new file mode 100644 index 00000000..a1ea1537 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C @@ -0,0 +1,772 @@ +/*---------------------------------------------------------------------------*\ +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(); + + tmp tnu = iter()*iter().nu(); + volScalarField& nu = tnu.ref(); + + for (++iter; iter != phases_.end(); ++iter) + { + nu += iter()*iter().nu(); + } + + 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 +{ + 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; +} + +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/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.H b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.H new file mode 100644 index 00000000..9a34c4ee --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.H @@ -0,0 +1,284 @@ +/*---------------------------------------------------------------------------*\ +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; + + 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/cfdemSolverMultiphase/multiphaseMixture/phase/phase.C b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.C new file mode 100644 index 00000000..ae68ff16 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.C @@ -0,0 +1,98 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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_) +{} + + +// * * * * * * * * * * * * * * * 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; + + if (nuModel_->read(phaseDict_)) + { + phaseDict_.lookup("rho") >> rho_; + + return true; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.H b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.H new file mode 100644 index 00000000..91341fcf --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.H @@ -0,0 +1,163 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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_; + + +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_; + } + + //- 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/cfdemSolverMultiphase/pEqn.H b/applications/solvers/cfdemSolverMultiphase/pEqn.H new file mode 100644 index 00000000..e5a374f0 --- /dev/null +++ b/applications/solvers/cfdemSolverMultiphase/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/UEqn.H b/applications/solvers/cfdemSolverPiso/UEqn.H index bfc3121d..2f26ac91 100644 --- a/applications/solvers/cfdemSolverPiso/UEqn.H +++ b/applications/solvers/cfdemSolverPiso/UEqn.H @@ -15,10 +15,10 @@ fvOptions.constrain(UEqn); if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull")) { solve(UEqn == - fvc::grad(p) + Ksl/rho*Us); - fvOptions.correct(U); + fvOptions.correct(U); } else if (piso.momentumPredictor()) { solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us); fvOptions.correct(U); -} \ No newline at end of file +} diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C index 8f1a8534..78c425c7 100644 --- a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C +++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C @@ -86,12 +86,12 @@ int main(int argc, char *argv[]) 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; + //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"); diff --git a/applications/solvers/cfdemSolverPiso/createFields.H b/applications/solvers/cfdemSolverPiso/createFields.H index 8f0b5584..36246301 100644 --- a/applications/solvers/cfdemSolverPiso/createFields.H +++ b/applications/solvers/cfdemSolverPiso/createFields.H @@ -96,17 +96,17 @@ #define createPhi_H Info<< "Reading/calculating face flux field phi\n" << endl; surfaceScalarField phi - ( - IOobject - ( +( + IOobject + ( "phi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE - ), - linearInterpolate(U*voidfraction) & mesh.Sf() - ); + ), + linearInterpolate(U*voidfraction) & mesh.Sf() +); #endif @@ -123,4 +123,4 @@ surfaceScalarField phi incompressible::turbulenceModel::New(U, phi, laminarTransport) ); -#include "createMRF.H" \ No newline at end of file +#include "createMRF.H" diff --git a/applications/solvers/cfdemSolverPiso/pEqn.H b/applications/solvers/cfdemSolverPiso/pEqn.H index 96ca1d77..c50eb9e3 100644 --- a/applications/solvers/cfdemSolverPiso/pEqn.H +++ b/applications/solvers/cfdemSolverPiso/pEqn.H @@ -31,12 +31,12 @@ constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF); while (piso.correctNonOrthogonal()) { // Pressure corrector - + fvScalarMatrix pEqn ( fvm::laplacian(rAUvoidfraction, p) == fvc::div(phi) + particleCloud.ddtVoidfraction() ); - + pEqn.setReference(pRefCell, pRefValue); pEqn.solve(mesh.solver(p.select(piso.finalInnerIter()))); @@ -55,4 +55,4 @@ else U = HbyA - voidfraction*rAU*fvc::grad(p) + Ksl/rho*Us*rAU; U.correctBoundaryConditions(); -fvOptions.correct(U); \ No newline at end of file +fvOptions.correct(U); diff --git a/applications/solvers/cfdemSolverPisoScalar/TEqn.H b/applications/solvers/cfdemSolverPisoScalar/TEqn.H index 6513b654..772b7838 100644 --- a/applications/solvers/cfdemSolverPisoScalar/TEqn.H +++ b/applications/solvers/cfdemSolverPisoScalar/TEqn.H @@ -1,4 +1,4 @@ - // get scalar source from DEM + // get scalar source from DEM particleCloud.forceM(1).manipulateScalarField(Tsource); Tsource.correctBoundaryConditions(); @@ -12,4 +12,4 @@ Tsource ); TEqn.relax(); - TEqn.solve(); \ No newline at end of file + TEqn.solve(); diff --git a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C index 9fc59b59..5f511cdf 100644 --- a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C +++ b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C @@ -81,23 +81,23 @@ int main(int argc, char *argv[]) { 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; + //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"); - + #include "TEqn.H" if(particleCloud.solveFlow()) diff --git a/applications/solvers/cfdemSolverPisoScalar/createFields.H b/applications/solvers/cfdemSolverPisoScalar/createFields.H index fe5dd29a..f7779bec 100644 --- a/applications/solvers/cfdemSolverPisoScalar/createFields.H +++ b/applications/solvers/cfdemSolverPisoScalar/createFields.H @@ -146,17 +146,17 @@ #define createPhi_H Info<< "Reading/calculating face flux field phi\n" << endl; surfaceScalarField phi - ( - IOobject - ( +( + IOobject + ( "phi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE - ), - linearInterpolate(U*voidfraction) & mesh.Sf() - ); + ), + linearInterpolate(U*voidfraction) & mesh.Sf() +); #endif @@ -173,4 +173,4 @@ surfaceScalarField phi incompressible::turbulenceModel::New(U, phi, laminarTransport) ); -#include "createMRF.H" \ No newline at end of file +#include "createMRF.H" diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index b2a30416..08af6a05 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H @@ -23,6 +23,9 @@ Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); + // correct source for the thermodynamic reference temperature + dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL)); + Qsource += QCoeff*Tref; fvScalarMatrix EEqn ( diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C index 39fa85cc..eb4d2c47 100644 --- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C +++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C @@ -69,8 +69,8 @@ int main(int argc, char *argv[]) #include "checkModelType.H" turbulence->validate(); - // #include "compressibleCourantNo.H" - // #include "setInitialDeltaT.H" + //#include "compressibleCourantNo.H" + //#include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/cfdemSolverRhoPimple/pEqn.H b/applications/solvers/cfdemSolverRhoPimple/pEqn.H index 40c68619..e107c5a4 100644 --- a/applications/solvers/cfdemSolverRhoPimple/pEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/pEqn.H @@ -32,7 +32,7 @@ else // + rhorAUf*fvc::ddtCorr(rho, U, phi) ) ); - + // flux without pressure gradient contribution phi = phiHbyA + phiUs; diff --git a/applications/solvers/cfdemSolverRhoPimple/rhoEqn.H b/applications/solvers/cfdemSolverRhoPimple/rhoEqn.H index 3be6fb2b..620e16c2 100644 --- a/applications/solvers/cfdemSolverRhoPimple/rhoEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/rhoEqn.H @@ -14,4 +14,4 @@ fvOptions.correct(rho); } -// ************************************************************************* // \ No newline at end of file +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H index edb4081f..5dfd1a5a 100644 --- a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H @@ -9,6 +9,10 @@ particleCloud.energyCoefficients(QCoeff); thCond=particleCloud.thermCondM().thermCond(); Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); +// correct source for the thermodynamic reference temperature +dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL)); +Qsource += QCoeff*Tref; + fvScalarMatrix EEqn ( fvm::ddt(rhoeps, he) + fvm::div(phi, he) diff --git a/applications/solvers/cfdemSolverRhoSimple/EEqn.H b/applications/solvers/cfdemSolverRhoSimple/EEqn.H index 3dbb5c3f..930f66e2 100644 --- a/applications/solvers/cfdemSolverRhoSimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/EEqn.H @@ -6,13 +6,13 @@ particleCloud.energyContributions(Qsource); particleCloud.energyCoefficients(QCoeff); - //thDiff=particleCloud.thermCondM().thermDiff(); - thCond=particleCloud.thermCondM().thermCond(); + //thDiff=particleCloud.thermCondM().thermDiff(); + thCond=particleCloud.thermCondM().thermCond(); - addSource = + addSource = ( he.name() == "e" - ? + ? fvc::div(phi, K) + fvc::div ( @@ -25,6 +25,9 @@ Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); + // correct source for the thermodynamic reference temperature + dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL)); + Qsource += QCoeff*Tref; fvScalarMatrix EEqn ( diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C index 4e1821fc..3af7f1f6 100644 --- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -17,14 +17,12 @@ License Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria Application - cfdemSolverRhoPimple + cfdemSolverRhoSimple Description - Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE) - algorithm. - Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. - The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x, - where additional functionality for CFD-DEM coupling is added. + Steady-state solver for turbulent flow of compressible fluids based on + rhoSimpleFoam where functionality for CFD-DEM coupling has been added. + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H index 5b2dbcbd..ed546344 100644 --- a/applications/solvers/cfdemSolverRhoSimple/pEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H @@ -27,7 +27,7 @@ else fvc::flux(rhoeps*HbyA) ) ); - + // flux without pressure gradient contribution phi = phiHbyA + phiUs; @@ -78,4 +78,4 @@ else U.correctBoundaryConditions(); fvOptions.correct(U); -K = 0.5*magSqr(U); \ No newline at end of file +K = 0.5*magSqr(U); diff --git a/applications/utilities/cfdemPostproc/cfdemPostproc.C b/applications/utilities/cfdemPostproc/cfdemPostproc.C index a832bacc..1f46b1de 100644 --- a/applications/utilities/cfdemPostproc/cfdemPostproc.C +++ b/applications/utilities/cfdemPostproc/cfdemPostproc.C @@ -84,7 +84,7 @@ int main(int argc, char *argv[]) particleCloud.dataExchangeM().allocateArray(particleV_,0.,1); particleCloud.get_cellIDs(cellIDs_); // get ref to cellIDs //particleCloud.dataExchangeM().allocateArray(cellIDs_,0.,1); - + while (runTime.loop()) { diff --git a/applications/utilities/writeUfluid/writeUfluid.C b/applications/utilities/writeUfluid/writeUfluid.C index cde6c88f..a6a84121 100644 --- a/applications/utilities/writeUfluid/writeUfluid.C +++ b/applications/utilities/writeUfluid/writeUfluid.C @@ -28,7 +28,7 @@ Application writeUfluidwriteUfluid Description - Writes the the cell center fluid velocity to particles in the lagrangian + Writes the the cell center fluid velocity to particles in the lagrangian time directory. \*---------------------------------------------------------------------------*/ @@ -76,13 +76,13 @@ int nParticle=0; { volVectorField U(UHeader,mesh); passiveParticleCloud myCloud(mesh, cloudName); - myCloud.write(); + myCloud.write(); nParticle = myCloud.size(); - IOField Ufluid(myCloud.fieldIOobject("Ufluid",IOobject::NO_READ),nParticle); + IOField Ufluid(myCloud.fieldIOobject("Ufluid",IOobject::NO_READ),nParticle); label i = 0; forAllConstIter(passiveParticleCloud, myCloud, iter) { - Ufluid[i]=U[iter().cell()]; + Ufluid[i]=U[iter().cell()]; i++; } Ufluid.write(); diff --git a/doc/CFDEMcoupling_install.txt b/doc/CFDEMcoupling_install.txt index 7878649b..b4634e3e 100644 --- a/doc/CFDEMcoupling_install.txt +++ b/doc/CFDEMcoupling_install.txt @@ -189,7 +189,7 @@ cfdemCompLIG :pre If the compilation fails with a message like -No rule to make target `/usr/lib/libpython2.7.so' :pre +No rule to make target '/usr/lib/libpython2.7.so' :pre you probably need to create a symbolic link to the library in question. diff --git a/doc/CFDEMcoupling_models.txt b/doc/CFDEMcoupling_models.txt index 246a4bb3..800e9eb4 100644 --- a/doc/CFDEMcoupling_models.txt +++ b/doc/CFDEMcoupling_models.txt @@ -38,7 +38,7 @@ models used for chemical reaction calculations. "diffusionCoefficients"_chemistryModel_diffusionCoefficients.html, "massTransferCoeff"_chemistryModel_massTransferCoeff.html, "off"_chemistryModel_noChemistry.html, -reactantPerParticle, +"reactantPerParticle"_chemistryModel_reactantPerParticle.html, "species"_chemistryModel_species.html :tb(c=2,ea=c) @@ -60,7 +60,8 @@ that performs the data exchange between the DEM code and the CFD code. "oneWayVTK"_dataExchangeModel_oneWayVTK.html, "twoWayFiles"_dataExchangeModel_twoWayFiles.html, "twoWayMPI"_dataExchangeModel_twoWayMPI.html, -"twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html :tb(c=2,ea=c) +"twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html, +"twoWayOne2One"_dataExchangeModel_twoWayOne2One.html :tb(c=2,ea=c) 6.6 Energy models :h4 @@ -94,11 +95,13 @@ Fines, "fieldStore"_forceModel_fieldStore.html, "fieldTimeAverage"_forceModel_fieldTimeAverage.html, "gradPForce"_forceModel_gradPForce.html, +"gradPForceSmooth"_forceModel_gradPForceSmooth.html, granKineticEnergy, "interface"_forceModel_interface.html, "noDrag"_forceModel_noDrag.html, "particleCellVolume"_forceModel_particleCellVolume.html, "pdCorrelation"_forceModel_pdCorrelation.html, +"surfaceTensionForce"_forceModel_surfaceTensionForce.html, "virtualMassForce"_forceModel_virtualMassForce.html, "viscForce"_forceModel_viscForce.html, "volWeightedAverage"_forceModel_volWeightedAverage.html :tb(c=2,ea=c) @@ -188,7 +191,8 @@ The "smoothingModel"_smoothingModel.html keyword entry specifies the model for smoothing the exchange fields. "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html, -"off"_smoothingModel_noSmoothing.html :tb(c=2,ea=c) +"off"_smoothingModel_noSmoothing.html, +"temporalSmoothing"_smoothingModel_temporalSmoothing.html :tb(c=2,ea=c) 6.16 Thermal conductivity models :h4 diff --git a/doc/CFDEMcoupling_solvers.txt b/doc/CFDEMcoupling_solvers.txt index 0f76f198..acbca05a 100644 --- a/doc/CFDEMcoupling_solvers.txt +++ b/doc/CFDEMcoupling_solvers.txt @@ -10,9 +10,10 @@ This section lists all CFDEMcoupling solvers alphabetically. "cfdemSolverIB"_cfdemSolverIB.html, +"cfdemSolverMultiphase"_cfdemSolverMultiphase.html, "cfdemSolverPiso"_cfdemSolverPiso.html, "cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html, -cfdemSolverRhoPimple, -cfdemSolverRhoPimpleChem, -cfdemSolverRhoSimple :tb(c=2,ea=c) +"cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html, +"cfdemSolverRhoPimpleChem"_cfdemSolverRhoPimpleChem.html, +"cfdemSolverRhoSimple"_cfdemSolverRhoSimple.html :tb(c=2,ea=c) diff --git a/doc/cfdemSolverMultiphase.txt b/doc/cfdemSolverMultiphase.txt new file mode 100644 index 00000000..6f7a08a3 --- /dev/null +++ b/doc/cfdemSolverMultiphase.txt @@ -0,0 +1,62 @@ + + + + + +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +cfdemSolverMultiphase command :h3 + +[Description:] + + +"cfdemSolverMultiphase" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam®(*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles. + + + + + +For more details, see "Vångö et al. (2018)"_#Vångö2018. + +:line + +:link(Vångö2018) +[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018): +"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds", +Applied Mathematical Modelling + +:line + + +NOTE: +(*) This offering is not approved or endorsed by OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com, and owner of the +OPENFOAM® and OpenCFD® trade marks. +OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com. + + + diff --git a/doc/cfdemSolverRhoPimple.txt b/doc/cfdemSolverRhoPimple.txt new file mode 100644 index 00000000..432581aa --- /dev/null +++ b/doc/cfdemSolverRhoPimple.txt @@ -0,0 +1,59 @@ + + + + + +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +cfdemSolverRhoPimple command :h3 + +[Description:] + + +"cfdemSolverRhoPimple" is a coupled CFD-DEM solver using the CFDEMcoupling +framework. Based on the OpenFOAM®(*) solver rhoPimpleFoam, this is a +transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE) +algorithm, coupled with the DEM code LIGGGHTS for solid particles. + + + + + +:line + + +NOTE: +(*) This offering is not approved or endorsed by OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com, and owner of the +OPENFOAM® and OpenCFD® trade marks. +OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com. + + + diff --git a/doc/cfdemSolverRhoPimpleChem.txt b/doc/cfdemSolverRhoPimpleChem.txt new file mode 100644 index 00000000..d3f64ddb --- /dev/null +++ b/doc/cfdemSolverRhoPimpleChem.txt @@ -0,0 +1,63 @@ + + + + + +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +cfdemSolverRhoPimpleChem command :h3 + +[Description:] + + +"cfdemSolverRhoPimpleChem" is a coupled CFD-DEM solver using the CFDEMcoupling +framework. Based on the OpenFOAM®(*) solver rhoPimpleFoam, this is a +transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE) +algorithm, coupled with the DEM code LIGGGHTS for solid particles. +Compared to cfdemSolverRhoPimple this solver adds functionality for chemical +reactions. + + + + + +:line + + +NOTE: +(*) This offering is not approved or endorsed by OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com, and owner of the +OPENFOAM® and OpenCFD® trade marks. +OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com. + + + diff --git a/doc/cfdemSolverRhoSimple.txt b/doc/cfdemSolverRhoSimple.txt new file mode 100644 index 00000000..d0ac3e00 --- /dev/null +++ b/doc/cfdemSolverRhoSimple.txt @@ -0,0 +1,59 @@ + + + + + +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +cfdemSolverRhoSimple command :h3 + +[Description:] + + +"cfdemSolverRhoSimple" is a coupled CFD-DEM solver using the CFDEMcoupling +framework. Based on the OpenFOAM®(*) solver rhoSimpleFoam, this is a +steady-state solver for turbulent flow of compressible fluids coupled with the +DEM code LIGGGHTS for solid particles. + + + + + +:line + + +NOTE: +(*) This offering is not approved or endorsed by OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com, and owner of the +OPENFOAM® and OpenCFD® trade marks. +OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com. + + + diff --git a/doc/chemistryModel_diffusionCoefficients.txt b/doc/chemistryModel_diffusionCoefficients.txt index df64243b..125a5584 100644 --- a/doc/chemistryModel_diffusionCoefficients.txt +++ b/doc/chemistryModel_diffusionCoefficients.txt @@ -23,7 +23,7 @@ diffusionCoefficientsProps diffusantGasNames ( speciesNames ); \} :pre -{switch1} = (optional, normally off) flag to give information :ulb,l +{switch1} = (optional, default false) flag to output verbose information :ulb,l {ChemistryFile} = path to file, where the reacting species are listed :l {diffusantGasNames} = list of gas field names that are the reactant gases :l :ule diff --git a/doc/chemistryModel_massTransferCoeff.txt b/doc/chemistryModel_massTransferCoeff.txt index 39573ac4..56534c52 100644 --- a/doc/chemistryModel_massTransferCoeff.txt +++ b/doc/chemistryModel_massTransferCoeff.txt @@ -21,7 +21,7 @@ massTransferCoeffProps verbose switch1; \} :pre -{switch1} = (optional, normally off) flag to give information :l +{switch1} = (optional, default false) flag to output verbose information :l :ule [Examples:] diff --git a/doc/chemistryModel_reactantPerParticle.txt b/doc/chemistryModel_reactantPerParticle.txt new file mode 100644 index 00000000..6f5069c3 --- /dev/null +++ b/doc/chemistryModel_reactantPerParticle.txt @@ -0,0 +1,54 @@ +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +chemistryModel reactantPerParticle command :h3 + +[Syntax:] + +Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties +dictionary. + +chemistryModels +( + reactantPerParticle +); +reactantPerParticleProps +\{ + voidfractionFieldName "voidfraction"; + Nevery number1; +\} :pre + +{voidfraction} = (optional, default "voidfraction") name of the finite volume void fraction field :l +{number1} = (optional, default 1) number to adjust execution interval :l +:ule + +[Examples:] + +chemistryModels +( + reactantPerParticle +); +reactantPerParticleProps +\{ + voidfractionFieldName "voidfraction"; + Nevery 1; +\} :pre + +[Description:] + +The chemistry model performs the calculation of chemical reactional effects +acting on each DEM particle. The reactantPerParticle model is the model to +communicate the available reactant per particle. + +[Restrictions:] + +none + +[Related commands:] + +"chemistryModel"_chemistryModel.html + diff --git a/doc/chemistryModel_species.txt b/doc/chemistryModel_species.txt index ad4891bf..b5071e59 100644 --- a/doc/chemistryModel_species.txt +++ b/doc/chemistryModel_species.txt @@ -26,16 +26,18 @@ speciesProps partTempName "partTemp"; partRhoName "partRho"; verbose switch1; + Nevery number1; \} :pre {ChemistryFile} = path to file, where the reacting species are listed :ulb,l -{T} = name of the finite volume temperature field, it is already added in default and doesn't need to be specified if name is the same :l -{rho} = name of the finite volume density field, it is already added in default and doesn't need to be specified if name is the same :l -{voidfraction} = name of the finite volume void fraction field, it is already added in default and doesn't need to be specified if name is the same :l -{molarConc} = name of the finite volume molar concentration field, it is already added in default and doesn't need to be specified if name is the same :l -{partTemp} = name of the finite volume cell averaged particle temperature field, it is already added in default and doesn't need to be specified if name is the same :l -{partRho} = name of the finite volume cell averaged density temperature field, it is already added in default and doesn't need to be specified if name is the same :l -{switch1} = (optional, normally off) flag to give information :l +{T} = (optional, default "T") name of the finite volume temperature field :l +{rho} = (optional, default "rho") name of the finite volume density field :l +{voidfraction} = (optional, default "voidfraction") name of the finite volume void fraction field :l +{molarConc} = (optional, default "molarConc") name of the finite volume molar concentration field :l +{partTemp} = (optional, default "partTemp") name of the finite volume cell averaged particle temperature field :l +{partRho} = (optional, default "partRho") name of the finite volume cell averaged density temperature field :l +{switch1} = (optional, default false) flag to output verbose information :l +{number1} = (optional, default 1) number to adjust execution interval :l :ule [Examples:] diff --git a/doc/dataExchangeModel_twoWayOne2One.txt b/doc/dataExchangeModel_twoWayOne2One.txt new file mode 100644 index 00000000..f415ebd5 --- /dev/null +++ b/doc/dataExchangeModel_twoWayOne2One.txt @@ -0,0 +1,47 @@ +"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c + +:link(lws,http://www.cfdem.com) +:link(lc,CFDEMcoupling_Manual.html#comm) + +:line + +dataExchangeModel_twoWayOne2One command :h3 + +[Syntax:] + +Defined in couplingProperties dictionary. + +dataExchangeModel twoWayOne2One; +twoWayOne2OneProps +\{ + liggghtsPath "path"; + useStaticProcMap switch1; +\}; :pre + +{path} = path to the DEM simulation input file :ulb,l +{switch1} = (optional, default no) switch to determine if the map is built once (yes) or every coupling step (no) :l +:ule + +[Examples:] + +dataExchangeModel twoWayOne2One; +twoWayOne2OneProps +\{ + liggghtsPath "../DEM/in.liggghts_init"; + useStaticProcMap yes; +\} :pre + +[Description:] + +The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayOne2One model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using a more sophisticated mapping scheme than twoWayMPI / all2all and scales much better for large systems and many cores. The DEM run is executed by the coupling model, via a liggghtsCommandModel object. Only use staticProcMap yes if no load balancing is employed. + +[Restrictions:] + +Must be used in combination with the engineSearchMany2Many locate model! Use the "one2one" cfd datacoupling option in fix couple/cfd in LIGGGHTS! + +Some warnings may be given for particles that have not been located - this is due to LIGGGHTS' treatment of domain crossers. + +[Related commands:] + +"dataExchangeModel"_dataExchangeModel.html + diff --git a/doc/forceModel_gradPForceSmooth.txt b/doc/forceModel_gradPForceSmooth.txt new file mode 100644 index 00000000..2b97c1ec --- /dev/null +++ b/doc/forceModel_gradPForceSmooth.txt @@ -0,0 +1,77 @@ +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +forceModel gradPForceSmooth command :h3 + +[Syntax:] + +Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties +dictionary. + +forceModels +( + gradPForceSmooth; +); +gradPForceSmoothProps +\{ + pFieldName "pressure"; + velocityFieldName "U"; + useAddedMass scalar1; + treatForceExplicit switch1; + treatForceDEM switch2; + interpolation switch3; + smoothingModel "smoothingModel"; +\} :pre + +{pressure} = name of the finite volume fluid pressure field :ulb,l +{U} = name of the finite volume fluid velocity field :l +{scalar1} = (optional, default 0) coefficient of added mass accounted for :l +{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l +{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l +{switch3} = (optional, default false) flag to use interpolated pressure values :l +{smoothingModel} = name of smoothing model :l +:ule + +[Examples:] + +forceModels +( + gradPForceSmooth; +); +gradPForceSmoothProps +\{ + pFieldName "p_rgh"; + velocityFieldName "U"; + interpolation false; + smoothingModel "temporalSmoothing"; + temporalSmoothingProps + \{ + lowerLimit 0.1; + upperLimit 1e10; + refField "p_rgh"; + gamma 1.0; + \} :pre + +\} :pre +[Description:] + +The {gradPForceSmooth} model calculates the particle based pressure gradient +force identically to the "gradPForce"_forceModel_gradPForce.html model but +allows smoothing of the pressure prior to the force calculation (without +altering the original pressure field). +Any smoothing model can be used and does not have to be the same as specified in +couplingProperties. Properties for the smoothing model have to be specified in a +sub-dictionary within {gradPForceSmoothProps}. + +[Restrictions:] + +A volScalarField "pSmooth" MUST be specified in the initial time directory! + +[Related commands:] + +"forceModel"_forceModel.html, "gradPForce"_forceModel_gradPForce.html + diff --git a/doc/forceModel_surfaceTensionForce.txt b/doc/forceModel_surfaceTensionForce.txt new file mode 100644 index 00000000..63657974 --- /dev/null +++ b/doc/forceModel_surfaceTensionForce.txt @@ -0,0 +1,55 @@ +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +forceModel surfaceTensionForce command :h3 + +[Syntax:] + +Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties +dictionary. + +forceModels +( + surfaceTensionForce; +); +surfaceTensionForceProps +\{ + stfFieldName "surfaceTensionField"; +\} :pre + +{surfaceTensionField} = name of the surface tension force field :ulb,l +:ule + +[Examples:] + +forceModels +( + surfaceTensionForce; +); +surfaceTensionForceProps +\{ + stfFieldName "surfaceTensionForce"; +\} :pre + +[Description:] + +The force model calculates the surface tension force acting on each DEM particle +based on a pre-calculated surface tension force as V_particle * F^sigma. When +used in conjunction with the "cfdemSolverMultiphase"_cfdemSolverMultiphase.html +solver, the surface tension force is calculated with the CSF (continuum surface +force) model (see Brackbill et al. (1992): "A continuum method for modeling +surface tension", J. Comput. Phys.). + +[Restrictions:] + +Has to be used with a multiphase solver that calculates the surface tension +force, e.g. {cfdemSolverMultiphase}. + +[Related commands:] + +"forceModel"_forceModel.html, "cfdemSolverMultiphase"_cfdemSolverMultiphase.html + diff --git a/doc/smoothingModel_temporalSmoothing.txt b/doc/smoothingModel_temporalSmoothing.txt new file mode 100644 index 00000000..eb6136e5 --- /dev/null +++ b/doc/smoothingModel_temporalSmoothing.txt @@ -0,0 +1,62 @@ +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +smoothingModel temporalSmoothing command :h3 + +[Syntax:] + +Defined in dictionary depending on the application. + +smoothingModel temporalSmoothing; +temporalSmoothingProps +\{ + lowerLimit number1; + upperLimit number2; + refField referenceField; + gamma smoothingStrength; +\} :pre + +{number1} = scalar fields will be bound to this lower value :ulb,l +{number2} = scalar fields will be bound to this upper value :l +{referenceField} = reference to the un-smoothed field required for the relaxation operation :l +{smoothingStrength} = control parameter for the smoothing, lower value yields stronger smoothing (gamma = 1 results in an equal contribution from the un-smoothed and smoothed fields) :l +:ule + +[Examples:] + +temporalSmoothingProps +\{ + lowerLimit 0.1; + upperLimit 1e10; + referenceField "p"; + gamma 1.0; +\} :pre + +[Description:] + +The {temporalSmoothing} model is a smoothing model that utilizes temporal +relaxation of a desired quantity. This model can be used to filter out +high-frequency fluctuations (e.g. numerical noise) controlled via the control +parameter gamma. +Note that this model does NOT smooth the calculated fields, instead smoothing is +performed on a separate (smooth) field which uses the calculated (un-smooth) +field as a reference. +Thus its usage is limited and CANNOT be used to smooth the exchange fields +similar to other smoothing models. +For further information see Vångö et al., "Unresolved CFD-DEM modeling of +multiphase flow in densely packed particle beds", Appl. Math. Model. (2018). + +[Restrictions:] + +This model does NOT smooth the calculated fields and can therefore NOT be used +as a general smoothing model to smoothen the exchange fields. +Attempting this will generate an error. + +[Related commands:] + +"smoothingModel"_smoothingModel.html + diff --git a/etc/bashrc b/etc/bashrc index e069a93e..c7224a0b 100755 --- a/etc/bashrc +++ b/etc/bashrc @@ -17,7 +17,7 @@ #------------------------------------------------------------------------------ export CFDEM_PROJECT=CFDEM -export CFDEM_VERSION=18.10 +export CFDEM_VERSION=19.02 ################################################################################ # USER EDITABLE PART: Changes made here may be lost with the next upgrade diff --git a/etc/cshrc b/etc/cshrc index b56a1c47..5ce854ac 100755 --- a/etc/cshrc +++ b/etc/cshrc @@ -15,7 +15,7 @@ #------------------------------------------------------------------------------ setenv CFDEM_PROJECT CFDEM -setenv CFDEM_VERSION 18.10 +setenv CFDEM_VERSION 19.02 ################################################################################ # USER EDITABLE PART: Changes made here may be lost with the next upgrade diff --git a/etc/library-list.txt b/etc/library-list.txt index e00d6b96..8e14f45d 100644 --- a/etc/library-list.txt +++ b/etc/library-list.txt @@ -1,3 +1,4 @@ lagrangian/cfdemParticle/dir lagrangian/cfdemParticleComp/dir finiteVolume/dir +../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 4bf3e896..c370a4d4 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -5,3 +5,4 @@ cfdemSolverRhoSimple/dir cfdemSolverIB/dir cfdemSolverPisoScalar/dir cfdemSolverRhoPimpleChem/dir +cfdemSolverMultiphase/dir diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index a184d5c1..2d409055 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -80,6 +80,8 @@ $(forceModels)/Fines/FanningDynFines.C $(forceModels)/Fines/ErgunStatFines.C $(forceModels)/granKineticEnergy/granKineticEnergy.C $(forceModels)/pdCorrelation/pdCorrelation.C +$(forceModels)/surfaceTensionForce/surfaceTensionForce.C +$(forceModels)/gradPForceSmooth/gradPForceSmooth.C $(forceModelsMS)/forceModelMS/forceModelMS.C $(forceModelsMS)/forceModelMS/newForceModelMS.C @@ -147,6 +149,8 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C $(dataExchangeModels)/twoWayFiles/twoWayFiles.C $(dataExchangeModels)/noDataExchange/noDataExchange.C $(dataExchangeModels)/twoWayMPI/twoWayMPI.C +$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C +$(dataExchangeModels)/twoWayOne2One/one2one.C $(averagingModels)/averagingModel/averagingModel.C $(averagingModels)/averagingModel/newAveragingModel.C @@ -169,5 +173,6 @@ $(smoothingModels)/smoothingModel/smoothingModel.C $(smoothingModels)/smoothingModel/newSmoothingModel.C $(smoothingModels)/noSmoothing/noSmoothing.C $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C +$(smoothingModels)/temporalSmoothing/temporalSmoothing.C LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME) diff --git a/src/lagrangian/cfdemParticle/Make/options b/src/lagrangian/cfdemParticle/Make/options index 9c152540..8a84bb35 100644 --- a/src/lagrangian/cfdemParticle/Make/options +++ b/src/lagrangian/cfdemParticle/Make/options @@ -34,6 +34,4 @@ LIB_LIBS = \ -lmpi_cxx \ -Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \ -L$(CFDEM_LIGGGHTS_BIN_DIR) \ - -lliggghts \ - -L$(CFDEM_Many2ManyLIB_PATH) \ - -lcoupleMany2Many + -lliggghts diff --git a/src/lagrangian/cfdemParticle/cfdTools/checkModelType.H b/src/lagrangian/cfdemParticle/cfdTools/checkModelType.H index a91bb786..a758d183 100644 --- a/src/lagrangian/cfdemParticle/cfdTools/checkModelType.H +++ b/src/lagrangian/cfdemParticle/cfdTools/checkModelType.H @@ -21,7 +21,7 @@ found=false; forAll(particleCloud.forceModels(),i) { - if(particleCloud.forceModels()[i]=="gradPForce") + if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth") found=true; } if(!found) @@ -56,7 +56,7 @@ found=false; forAll(particleCloud.forceModels(),i) { - if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="viscForce") + if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth" || particleCloud.forceModels()[i]=="viscForce") found=true; } if(found) @@ -80,7 +80,7 @@ found=false; forAll(particleCloud.forceModels(),i) { - if(particleCloud.forceModels()[i]=="gradPForce") + if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth") found=true; } if(!found) @@ -99,3 +99,6 @@ Warning << "You chose model type -none- you might get erroneous results!" << endl; else FatalError <<"no suitable model type specified:" << modelType << "\n" << abort(FatalError); + + if (particleCloud.smoothingM().type() == "temporalSmoothing") + FatalError << "the temporalSmoothing model does not support smoothing of the exchange fields, please see documentation!" << endl; diff --git a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H index 637262b3..7297ac62 100755 --- a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H +++ b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H @@ -34,8 +34,8 @@ Description #ifndef versionInfo_H #define versionInfo_H -word CFDEMversion="PFM 18.10"; -word compatibleLIGGGHTSversion="PFM 18.10"; +word CFDEMversion="PFM 19.02"; +word compatibleLIGGGHTSversion="PFM 19.02"; word OFversion="4.x"; Info << "\nCFDEMcoupling version: " << CFDEMversion << endl; diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C index f9a2eb9d..62695fd1 100644 --- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C +++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C @@ -312,7 +312,7 @@ void diffusionCoefficient::execute() TotalFraction_[i] += Xfluid_[j]/dBinary_; // dCoeff -- diffusion component of diffusant gas - MixtureBinaryDiffusion_[i] = (1.0-XfluidDiffusant_[i])/TotalFraction_[i]; + MixtureBinaryDiffusion_[i] = (1.0-XfluidDiffusant_[i])/TotalFraction_[i]; if(verbose_) { diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H index fb96c5f8..ab360a67 100644 --- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H +++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H @@ -72,7 +72,7 @@ private: // gas pressure at particle location word pressureFieldName_; - const volScalarField& P_; + const volScalarField& P_; word partPressureName_; diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C index 57a58f44..1573c68a 100644 --- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C +++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C @@ -109,8 +109,6 @@ void reactantPerParticle::reAllocMyArrays() const void reactantPerParticle::execute() { - - loopCounter_++; if (loopCounter_ % Nevery_ != 0) { @@ -121,7 +119,7 @@ void reactantPerParticle::execute() particlesPerCell_ *= 0.0; - label cellI=0; + label cellI=0; scalar voidfraction(1); scalar cellvolume(0.0); scalar particlesPerCell(1.0); @@ -133,7 +131,7 @@ void reactantPerParticle::execute() if (cellI >= 0) { particlesPerCell_[cellI] += 1.0; - } + } } // no fill array and communicate it @@ -151,7 +149,7 @@ void reactantPerParticle::execute() // give DEM data particleCloud_.dataExchangeM().giveData("reactantPerParticle", "scalar-atom", reactantPerParticle_); - + Info << "give data done" << endl; } diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H index 6153b593..672f1138 100755 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H @@ -245,17 +245,25 @@ public: virtual int getNumberOfTypes() const; virtual double* getTypeVol() const; - inline void setPositions(label n,double* pos) const + inline void setPositions(label n,double* pos) { for (int i=0;i. + +Copyright 2018- Paul Kieckhefen, TUHH + +\*---------------------------------------------------------------------------*/ + +//#define O2ODEBUG + +#ifdef O2ODEBUG +#include +#endif + +#include +#include "one2one.H" + +template +struct mpi_type_wrapper { + MPI_Datatype mpi_type; + mpi_type_wrapper(); +}; +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_FLOAT) {} +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_DOUBLE) {} +template <> mpi_type_wrapper::mpi_type_wrapper() +: mpi_type(MPI_INT) {} + +/* ---------------------------------------------------------------------- */ + +One2One::One2One(MPI_Comm caller) +: +ncollected_(-1), +comm_(caller), +nsrc_procs_(-1), +src_procs_(nullptr), +ndst_procs_(-1), +dst_procs_(nullptr), +nlocal_(-1), +natoms_(nullptr), +request_(nullptr), +status_(nullptr) +{ + MPI_Comm_rank(comm_,&me_); + MPI_Comm_size(comm_,&nprocs_); +} + +/* ---------------------------------------------------------------------- */ + +One2One::~One2One() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- + communicate particle ids based on processor communication pattern +------------------------------------------------------------------------- */ + +void One2One::setup +( + int nsrc_procs, + int *src_procs, + int ndst_procs, + int *dst_procs, + int nlocal +) +{ + // free any previous one2one info + deallocate(); + + src_procs_ = src_procs; + nsrc_procs_ = nsrc_procs; + dst_procs_ = dst_procs; + ndst_procs_ = ndst_procs; + nlocal_ = nlocal; + + // gather number of ids for reserving memory + natoms_ = new int[nprocs_]; + MPI_Allgather // may be replaced by send/irecv + ( + &nlocal_, + 1, + MPI_INT, + natoms_, + 1, + MPI_INT, + comm_ + ); + + ncollected_ = 0; + int nrequests = 0; + for (int i = 0; i < nsrc_procs_; i++) + { + if (natoms_[src_procs_[i]] > 0) + { + ncollected_ += natoms_[src_procs_[i]]; + + if (src_procs_[i] != me_) // no receive for on-proc info + { + nrequests++; + } + } + } + + if (nrequests > 0) + { + request_ = new MPI_Request[nrequests]; + status_ = new MPI_Status[nrequests]; + } +} + +/* ---------------------------------------------------------------------- + src: what is present on this proc + dst: what is received from other procs + all comm according to map set up in setup(...) +------------------------------------------------------------------------- */ +template +void One2One::exchange(T *&src, T *&dst, int data_length) +{ + mpi_type_wrapper wrap; + + // post receives + int offset_local = -1; + int offset = 0; + int requesti = 0; + for (int i = 0; i < nsrc_procs_; i++) + { + // do post a receives for procs who own particles + if (natoms_[src_procs_[i]] > 0) + { + if (src_procs_[i] != me_) + { + #ifdef O2ODEBUG + std::cout<< "[" << me_ << "]" + << " RCV " << i + << " of " << nsrc_procs_ + << " from: " << src_procs_[i] + << " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]] + << " datalength " << data_length + << " offset " << offset + << std::endl; + #endif + MPI_Irecv + ( + &dst[offset], + natoms_[src_procs_[i]]*data_length, + wrap.mpi_type, + src_procs_[i], + MPI_ANY_TAG, + comm_, + &request_[requesti] + ); + requesti++; + } + else // data is available on-proc + { + offset_local = offset; + } + } + offset += natoms_[src_procs_[i]]*data_length; + } + + // make sure all receives are posted + MPI_Barrier(comm_); + + // blocking sends - do nonblocking instead + // since doing many-2-many here? + // only do sends if I have particles + if (nlocal_ > 0) + { + for (int i = 0; i < ndst_procs_; i++) + { + if (dst_procs_[i] != me_) + { + #ifdef O2ODEBUG + std::cout<< "[" << me_ << "]" + << " SEND to: " << dst_procs_[i] + << " nlocal_ " << nlocal_ + << " data_length " << data_length + << std::endl; + #endif + MPI_Send + ( + src, + nlocal_*data_length, + wrap.mpi_type, + dst_procs_[i], + 0, + comm_ + ); + } + } + } + + // only wait if requests were actually posted + if (requesti > 0) + MPI_Waitall(requesti, request_, status_); + + // copy on-proc data + if (offset_local > -1) + { + const int max_locali = nlocal_ * data_length; + for + ( + int locali = 0; + locali < max_locali; + locali++ + ) + { + dst[locali+offset_local] = src[locali]; + } + } +} + +template void One2One::exchange(int*&, int*&, int); +template void One2One::exchange(double*&, double*&, int); + +// there should be a way to do this without copying data +template +void One2One::exchange(T **&src, T **&dst, int data_length) +{ + mpi_type_wrapper wrap; + + T* tmp_dst = new T[ncollected_*data_length]; + T* tmp_src = new T[nlocal_*data_length]; + + for (int i = 0; i < nlocal_; i++) + for (int j = 0; j < data_length; j++) + tmp_src[data_length*i+j] = src[i][j]; + + exchange(tmp_src, tmp_dst, data_length); + + for (int i = 0; i < ncollected_; i++) + for (int j = 0; j < data_length; j++) + dst[i][j] = tmp_dst[data_length*i+j]; + + delete [] tmp_src; + delete [] tmp_dst; +} +template void One2One::exchange(int**&, int**&, int); +template void One2One::exchange(double**&, double**&, int); + + +template +void One2One::exchange(T **&src, T *&dst, int data_length) +{ + mpi_type_wrapper wrap; + + T* tmp_src = new T[nlocal_*data_length]; + + for (int i = 0; i < nlocal_; i++) + for (int j = 0; j < data_length; j++) + tmp_src[data_length*i+j] = src[i][j]; + + exchange(tmp_src, dst, data_length); + + + delete [] tmp_src; +} +template void One2One::exchange(int**&, int*&, int); +template void One2One::exchange(double**&, double*&, int); + +/* ---------------------------------------------------------------------- */ + +void One2One::deallocate() +{ + delete [] src_procs_; + delete [] dst_procs_; + delete [] natoms_; + + delete [] request_; + delete [] status_; +} diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H new file mode 100644 index 00000000..a3d1e21e --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ +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 2018- Paul Kieckhefen, TUHH + +\*---------------------------------------------------------------------------*/ + +#ifndef ONE2ONE_H +#define ONE2ONE_H + +#include + +class One2One +{ + public: + One2One(MPI_Comm); + ~One2One(); + + void setup + ( + int nsrc_procs, + int *src_procs, + int ndst_procs, + int* dst_procs, + int nlocal + ); + + template + void exchange(T *&, T *&, int data_length=1); + template + void exchange(T **&, T **&, int data_length=1); + template + void exchange(T **&, T *&, int data_length=1); + + int ncollected_; // # of ids in from group + + protected: + + MPI_Comm comm_; + int me_, nprocs_; // rank and size + + // communication partners + int nsrc_procs_; // # of off-processor IDs + int* src_procs_; // procs I receive data from + int ndst_procs_; // # of off-processor IDs + int* dst_procs_; // procs I receive data from + + int nlocal_; // # particle ids I own + int* natoms_; + + MPI_Request* request_; + MPI_Status* status_; + + void deallocate(); +}; + +#endif diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C new file mode 100644 index 00000000..d5738c59 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C @@ -0,0 +1,941 @@ +/*---------------------------------------------------------------------------*\ + 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). + +Contributing authors + Paul Kieckhefen (TUHH) 2018- + +\*---------------------------------------------------------------------------*/ + + +#include "twoWayOne2One.H" +#include "addToRunTimeSelectionTable.H" +#include "clockModel.H" +#include "pair.h" +#include "force.h" +#include "forceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(twoWayOne2One, 0); + +addToRunTimeSelectionTable +( + dataExchangeModel, + twoWayOne2One, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +twoWayOne2One::twoWayOne2One +( + const dictionary& dict, + cfdemCloud& sm +) +: + dataExchangeModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + thisLigPartner_(0), + thisFoamPartner_(0), + lig2foam_(nullptr), + foam2lig_(nullptr), + lig2foam_mask_(nullptr), + lig2foam_ids_(nullptr), + foam2lig_ids_(nullptr), + lig2foam_vec_tmp_(nullptr), + lig2foam_scl_tmp_(nullptr), + foam2lig_vec_tmp_(nullptr), + foam2lig_scl_tmp_(nullptr), + staticProcMap_(propsDict_.lookupOrDefault("useStaticProcMap", false)), + cellIdComm_(propsDict_.lookupOrDefault("useCellIdComm", false)), + my_prev_cell_ids_fix_(nullptr), + verbose_(propsDict_.lookupOrDefault("verbose", false)), + lmp(nullptr) +{ + Info<<"Starting up LIGGGHTS for first time execution"<input->file(liggghtsPath.c_str()); + + // get DEM time step size + DEMts_ = lmp->update->dt; + checkTSsize(); + + // calculate boundingBox of FOAM subdomain + primitivePatch tmpBoundaryFaces + ( + SubList + ( + sm.mesh().faces(), + sm.mesh().nFaces() - sm.mesh().nInternalFaces(), + sm.mesh().nInternalFaces() + ), + sm.mesh().points() + ); + typedef PrimitivePatch bPatch; + bPatch boundaryFaces + ( + tmpBoundaryFaces.localFaces(), + tmpBoundaryFaces.localPoints() + ); + thisFoamBox_ = treeBoundBox(boundaryFaces.localPoints()); + if (staticProcMap_) + { + createProcMap(); + } + + if (cellIdComm_) + { + my_prev_cell_ids_fix_ = static_cast + ( lmp->modify->find_fix_property + ( + "prev_cell_ids", + "property/atom", + "scalar", + 0, + 0, + "cfd coupling", + true + ) + ); + } +} + +void twoWayOne2One::createProcMap() +{ + List foamBoxes(Pstream::nProcs()); + foamBoxes[Pstream::myProcNo()] = thisFoamBox_; + Pstream::gatherList(foamBoxes); + Pstream::scatterList(foamBoxes); + + // calculate bounding box of LIG subdomain + // this may have to move to couple when dynamic LB occurs + List ligBoxes(Pstream::nProcs()); + double** ligbb = o2o_liggghts_get_boundingbox(lmp); + boundBox thisLigBox + ( + point(ligbb[0][0], ligbb[0][1], ligbb[0][2]), + point(ligbb[1][0], ligbb[1][1], ligbb[1][2]) + ); + ligBoxes[Pstream::myProcNo()] = thisLigBox; + Pstream::gatherList(ligBoxes); + Pstream::scatterList(ligBoxes); + + thisLigPartner_.clear(); + thisFoamPartner_.clear(); + + // detect LIG subdomains which this FOAM has to interact with + forAll(ligBoxes, ligproci) + { + if (thisFoamBox_.overlaps(ligBoxes[ligproci])) + { + thisLigPartner_.append(ligproci); + } + } + // detect FOAM subdomains this LIG has to interact with + // TODO: refactor to invert this list here + forAll(foamBoxes, foamproci) + { + if (thisLigBox.overlaps(foamBoxes[foamproci])) + { + thisFoamPartner_.append(foamproci); + } + } + + if (verbose_) + { + Pout<< "FOAM bounding box: " << thisFoamBox_ + << " LIG bounding box: " << thisLigBox + << nl + << "FOAM comm partners: " << thisFoamPartner_ + << " LIG comm partners: " << thisLigPartner_ + << endl; + } +} + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +twoWayOne2One::~twoWayOne2One() +{ + delete foam2lig_; + delete lig2foam_; + + destroy(lig2foam_ids_); + destroy(foam2lig_ids_); + + destroy(lig2foam_vec_tmp_); + destroy(lig2foam_scl_tmp_); + destroy(foam2lig_vec_tmp_); + destroy(foam2lig_scl_tmp_); + + delete lmp; +} + +// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // +void twoWayOne2One::getData +( + word name, + word type, + double ** const& field, + label /*step*/ +) const +{ + if (name == "x") // the location is transferred by couple() + { + return; + } + if (type == "vector-atom") + { + double **tmp_= static_cast(lammps_extract_atom(lmp,name.c_str())); + if (!tmp_) + { + LAMMPS_NS::Fix *fix = nullptr; + fix = lmp->modify->find_fix_property + ( + name.c_str(), + "property/atom", + "vector", + 0, + 0, + "cfd coupling", + false + ); + if (fix) + { + tmp_ = static_cast(fix)->array_atom; + } + else + { + Warning<< "coupling fix not found!" << endl; + } + + if (!tmp_) + { + FatalError<< "find_fix_property " << name + << " array_atom not found." + << abort(FatalError); + } + } + + lig2foam_->exchange + ( + tmp_, + lig2foam_vec_tmp_, + 3 + ); + extractCollected + ( + lig2foam_vec_tmp_, + const_cast(field), + 3 + ); + } + else if (type == "scalar-atom") + { + double *tmp_ = static_cast(lammps_extract_atom(lmp,name.c_str())); + if (!tmp_) + { + LAMMPS_NS::Fix *fix = nullptr; + fix = lmp->modify->find_fix_property + ( + name.c_str(), + "property/atom", + "scalar", + 0, + 0, + "cfd coupling", + true + ); + + if (fix) + { + tmp_ = static_cast(fix)->vector_atom; + } + else + { + FatalError<< "coupling fix not found!" << abort(FatalError); + } + + if (!tmp_) + { + FatalError<< "find_fix_property " << name + << " vector_atom not found." + << abort(FatalError); + } + } + lig2foam_->exchange + ( + tmp_, + lig2foam_scl_tmp_ + ); + extractCollected + ( + lig2foam_scl_tmp_, + const_cast(field) + ); + } + else + { + FatalError << "requesting type " << type << " and name " << name << abort(FatalError); + } +} + +void twoWayOne2One::getData +( + word name, + word type, + int ** const& field, + label /*step*/ +) const +{ + FatalError << "do not use this getData!!!" << abort(FatalError); +/* + o2o_data_liggghts_to_of + ( + name.c_str(), + type.c_str(), + lmp, + (void*&) field, + "int", + comm_liggghts_ + ); +*/ +} + +void twoWayOne2One::giveData +( + word name, + word type, + double ** const& field, + const char* datatype +) const +{ + if (type == "vector-atom") + { + foam2lig_->exchange + ( + const_cast(field), + foam2lig_vec_tmp_, + 3 + ); + o2o_data_of_to_liggghts + ( + name.c_str(), + type.c_str(), + lmp, + foam2lig_vec_tmp_, + datatype, + foam2lig_ids_, + foam2lig_->ncollected_ + ); + } + else if (type == "scalar-atom") + { + foam2lig_->exchange + ( + const_cast(field), + foam2lig_scl_tmp_, + 1 + ); + o2o_data_of_to_liggghts + ( + name.c_str(), + type.c_str(), + lmp, + foam2lig_scl_tmp_, + datatype, + foam2lig_ids_, + foam2lig_->ncollected_ + ); + } + else + { + FatalError<< "twoWayMany2Many::giveData requested type " << type + << " not implemented!" + << abort(FatalError); + } +} + +//============ +// double ** +void twoWayOne2One::allocateArray +( + double**& array, + double initVal, + int width, + int length +) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, width, "o2o:dbl**"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void twoWayOne2One::allocateArray +( + double**& array, + double initVal, + int width, + const char* length +) const +{ + int len = max(particleCloud_.numberOfParticles(),1); + lmp->memory->grow(array, len, width, "o2o:dbl**:autolen"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void inline twoWayOne2One::destroy(double** array,int len) const +{ + lmp->memory->destroy(array); +} + +//============ +// int ** +void twoWayOne2One::allocateArray +( + int**& array, + int initVal, + int width, + int length +) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, width, "o2o:int**"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void twoWayOne2One::allocateArray +( + int**& array, + int initVal, + int width, + const char* length +) const +{ + int len = max(particleCloud_.numberOfParticles(),1); + lmp->memory->grow(array, len, width, "o2o:int**:autolen"); + for (int i = 0; i < len; i++) + for (int j = 0; j < width; j++) + array[i][j] = initVal; +} + +void inline twoWayOne2One::destroy(int** array,int len) const +{ + lmp->memory->destroy(array); +} + +//============ +// double * +void twoWayOne2One::allocateArray(double*& array, double initVal, int length) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, "o2o:dbl*"); + for (int i = 0; i < len; i++) + array[i] = initVal; +} + +void inline twoWayOne2One::destroy(double* array) const +{ + lmp->memory->destroy(array); +} + +//============== +// int * +void twoWayOne2One::allocateArray(int*& array, int initVal, int length) const +{ + int len = max(length,1); + lmp->memory->grow(array, len, "o2o:int*"); + for (int i = 0; i < len; i++) + array[i] = initVal; +} + +void inline twoWayOne2One::destroy(int* array) const +{ + lmp->memory->destroy(array); +} +//============== + +bool twoWayOne2One::couple(int i) +{ + bool coupleNow = false; + if (i==0) + { + couplingStep_++; + coupleNow = true; + + + // run commands from liggghtsCommands dict + Info<< "Starting up LIGGGHTS" << endl; + particleCloud_.clockM().start(3,"LIGGGHTS"); + + // check if liggghtsCommandModels with exaxt timing are being run + bool exactTiming(false); + int runComNr = -10; + DynamicList interruptTimes(0); + DynamicList DEMstepsToInterrupt(0); + DynamicList lcModel(0); + + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + // Check if exact timing is needed + // get time for execution + // store time for execution in list + if(particleCloud_.liggghtsCommand(i).exactTiming()) + { + exactTiming = true; + DynamicList h + = particleCloud_.liggghtsCommand(i).executionsWithinPeriod + ( + TSstart(), + TSend() + ); + + forAll(h,j) + { + // save interrupt times (is this necessary) + interruptTimes.append(h[j]); + + // calc stepsToInterrupt + DEMstepsToInterrupt.append(DEMstepsTillT(h[j])); + + // remember which liggghtsCommandModel to run + lcModel.append(i); + } + + // make cumulative + label len = DEMstepsToInterrupt.size(); + label ind(0); + forAll(DEMstepsToInterrupt,i) + { + ind = len - i - 1; + if(ind > 0) + { + DEMstepsToInterrupt[ind] -= DEMstepsToInterrupt[ind-1]; + } + } + + Info<< "Foam::twoWayOne2One::couple(i): interruptTimes=" << interruptTimes << nl + << "Foam::twoWayOne2One::couple(i): DEMstepsToInterrupt=" << DEMstepsToInterrupt << nl + << "Foam::twoWayOne2One::couple(i): lcModel=" << lcModel + << endl; + } + + if(particleCloud_.liggghtsCommand(i).type() == "runLiggghts") + { + runComNr = i; + } + } + + // models with exact timing exists + label commandLines(0); + if(exactTiming) + { + // extension for more liggghtsCommands active the same time: + // sort interrupt list within this run period + // keep track of corresponding liggghtsCommand + int DEMstepsRun(0); + + forAll(interruptTimes,j) + { + // set run command till interrupt + DEMstepsRun += DEMstepsToInterrupt[j]; + particleCloud_.liggghtsCommand(runComNr).set(DEMstepsToInterrupt[j]); + const char* command = particleCloud_.liggghtsCommand(runComNr).command(0); + Info<< "Executing run command: '"<< command <<"'"<< endl; + lmp->input->one(command); + + // run liggghts command with exact timing + command = particleCloud_.liggghtsCommand(lcModel[j]).command(0); + Info << "Executing command: '"<< command <<"'"<< endl; + lmp->input->one(command); + } + + // do the run + if(particleCloud_.liggghtsCommand(runComNr).runCommand(couplingStep())) + { + particleCloud_.liggghtsCommand(runComNr).set(couplingInterval() - DEMstepsRun); + const char* command = particleCloud_.liggghtsCommand(runComNr).command(0); + Info<< "Executing run command: '"<< command <<"'"<< endl; + lmp->input->one(command); + } + + // do the other non exact timing models + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + if + ( + ! particleCloud_.liggghtsCommand(i).exactTiming() && + particleCloud_.liggghtsCommand(i).runCommand(couplingStep()) + ) + { + commandLines=particleCloud_.liggghtsCommand(i).commandLines(); + for(int j=0;jinput->one(command); + } + } + } + } + // no model with exact timing exists + else + { + forAll(particleCloud_.liggghtsCommandModelList(),i) + { + if(particleCloud_.liggghtsCommand(i).runCommand(couplingStep())) + { + commandLines=particleCloud_.liggghtsCommand(i).commandLines(); + for(int j=0;jinput->one(command); + } + } + } + } + + particleCloud_.clockM().stop("LIGGGHTS"); + Info<< "LIGGGHTS finished" << endl; + + if (!staticProcMap_) + { + createProcMap(); + } + + setupLig2FoamCommunication(); + + locateParticles(); + + setupFoam2LigCommunication(); + + if (verbose_) + { + Pout<< "FOAM owns " << getNumberOfParticles() + << " LIG owns " << lmp->atom->nlocal + << nl + << "FOAM collects " << lig2foam_->ncollected_ + << " LIG collects " << foam2lig_->ncollected_ + << endl; + } + } + + return coupleNow; +} + +void twoWayOne2One::setupLig2FoamCommunication() +{ + int* src_procs = new int[thisLigPartner_.size()]; + for (int proci = 0; proci < thisLigPartner_.size(); proci++) + { + src_procs[proci] = thisLigPartner_[proci]; + } + int* dst_procs = new int[thisFoamPartner_.size()]; + for (int proci = 0; proci < thisFoamPartner_.size(); proci++) + { + dst_procs[proci] = thisFoamPartner_[proci]; + } + + delete lig2foam_; + lig2foam_ = new One2One(comm_liggghts_); + lig2foam_->setup + ( + thisLigPartner_.size(), + src_procs, + thisFoamPartner_.size(), + dst_procs, + lmp->atom->nlocal + ); + allocateArray + ( + lig2foam_vec_tmp_, + 0., + 3 * lig2foam_->ncollected_ + ); + allocateArray + ( + lig2foam_scl_tmp_, + 0., + lig2foam_->ncollected_ + ); +} + + +void twoWayOne2One::locateParticles() +{ + // get positions for locate + double** my_positions = static_cast(lmp->atom->x); + double* my_flattened_positions = nullptr; + allocateArray(my_flattened_positions, 0., 3*lmp->atom->nlocal); + for (int atomi = 0; atomi < lmp->atom->nlocal; atomi++) + { + for (int coordi = 0; coordi < 3; coordi++) + { + my_flattened_positions[atomi*3+coordi] = my_positions[atomi][coordi]; + } + } + + double* collected_flattened_positions = nullptr; + allocateArray(collected_flattened_positions, 0., 3*lig2foam_->ncollected_); + + lig2foam_->exchange(my_flattened_positions, collected_flattened_positions, 3); + destroy(my_flattened_positions); + + double* my_prev_cell_ids = nullptr; + double* prev_cell_ids = nullptr; + if (cellIdComm_) + { + my_prev_cell_ids = my_prev_cell_ids_fix_->vector_atom; + allocateArray(prev_cell_ids, -1, lig2foam_->ncollected_); + lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids); + } + + if (lig2foam_mask_) + { + delete [] lig2foam_mask_; + } + lig2foam_mask_ = new bool[lig2foam_->ncollected_]; + + DynamicList