zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
mirror of https://github.com/ParticulateFlow/CFDEMcoupling-PFM.git synced 2025-12-08 06:37:44 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge pull request #146 from ParticulateFlow/release

Browse Source 
Release 23.02
This commit is contained in:
Daniel
2023-02-15 12:26:45 +01:00
committed by GitHub
parent 16a415a152 13e5989216
commit 86753f2823
316 changed files with 12811 additions and 22422 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -8,7 +8,7 @@ log.*
*.swp
*.swo
**/linux*Gcc*/
**/linux*cc*/
**/.vscode
lnInclude

3
applications/solvers/cfdemSolverIB/Make/options
View File

@ -9,12 +9,12 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \
@ -27,6 +27,7 @@ EXE_LIBS = \
-ldynamicFvMesh \
-ldynamicMesh \
-lfvOptions \
-lsampling \
-l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS)

33
applications/solvers/cfdemSolverIB/cfdemSolverIB.C
View File

@ -6,7 +6,8 @@
Christoph Goniva, christoph.goniva@cfdem.com
Copyright (C) 1991-2009 OpenCFD Ltd.
Copyright (C) 2009-2012 JKU, Linz
Copyright (C) 2012- DCS Computing GmbH,Linz
Copyright (C) 2012-2015 DCS Computing GmbH,Linz
Copyright (C) 2015- JKU, Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling.
@ -29,11 +30,14 @@ Application
Description
Transient solver for incompressible flow.
The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
where additional functionality for CFD-DEM coupling using immersed body
(fictitious domain) method is added.
Contributions
Alice Hager
Daniel Queteschiner
Thomas Lichtenegger
Achuth N. Balachandran Nair
\*---------------------------------------------------------------------------*/
@ -53,23 +57,21 @@ Contributions
#include "cellSet.H"
#include "fvOptions.H" // added the fvOptions library
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "createFvOptions.H"
// create cfdemCloud
#include "readGravitationalAcceleration.H"
@ -93,24 +95,31 @@ int main(int argc, char *argv[])
// do particle stuff
Info << "- evolve()" << endl;
particleCloud.evolve();
particleCloud.evolve(Us);
// Pressure-velocity PISO corrector
{
MRF.correctBoundaryVelocity(U);
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U)
fvm::ddt(voidfraction,U) + MRF.DDt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}
// --- PISO loop
@ -126,6 +135,7 @@ int main(int argc, char *argv[])
adjustPhi(phi, U, p);
while (piso.correctNonOrthogonal())
{
// Pressure corrector
@ -152,12 +162,15 @@ int main(int argc, char *argv[])
}
}
laminarTransport.correct();
turbulence->correct();
Info << "particleCloud.calcVelocityCorrection() " << endl;
volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
fvOptions.correct(U);
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

35
applications/solvers/cfdemSolverIB/createFields.H
View File

@ -26,21 +26,6 @@
),
mesh
);
//mod by alice
Info<< "Reading physical velocity field U" << endl;
Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// drag law modelling
@ -76,9 +61,8 @@
mesh
);
//mod by alice
Info<< "Reading field phiIB\n" << endl;
Info<< "Reading field voidfraction\n" << endl;
volScalarField voidfraction
(
IOobject
@ -91,6 +75,21 @@
),
mesh
);
Info<< "Reading particle velocity field Us\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
# include "createPhi.H"
@ -126,3 +125,5 @@
);
//===========================
#include "createMRF.H"

3
applications/solvers/cfdemSolverIBContinuousForcing/Make/files Executable file
View File

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

34
applications/solvers/cfdemSolverIBContinuousForcing/Make/options Executable file
View File

@ -0,0 +1,34 @@
include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
EXE_INC = \
-I$(CFDEM_OFVERSION_DIR) \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools \
-ldynamicFvMesh \
-ldynamicMesh \
-lfvOptions \
-lsampling \
-l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS)

19
applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H Normal file
View File

@ -0,0 +1,19 @@
fvVectorMatrix UEqn
(
fvm::ddt(voidfractionNext,U) + MRF.DDt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
+ (lambda*(1-voidfractionNext)/U.mesh().time().deltaT())*(fvc::Sp(1,Us)-fvm::Sp(1,U))
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}

129
applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C Executable file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling - Open Source CFD-DEM coupling
CFDEMcoupling is part of the CFDEMproject
www.cfdem.com
Copyright (C) 1991-2009 OpenCFD Ltd.
Copyright (C) 2009-2012 JKU, Linz
Copyright (C) 2012-2015 DCS Computing GmbH,Linz
Copyright (C) 2015- JKU, Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling.
CFDEMcoupling is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with CFDEMcoupling. If not, see <http://www.gnu.org/licenses/>.
Application
cfdemSolverIBContinuousForcing
Description
Transient solver for incompressible flow.
The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
where additional functionality for CFD-DEM coupling using immersed body
(fictitious domain) method and a continuous forcing approach is added.
Contributions
Alice Hager
Achuth N. Balachandran Nair
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "cfdemCloudIB.H"
#include "implicitCouple.H"
#include "averagingModel.H"
#include "regionModel.H"
#include "voidFractionModel.H"
#include "dynamicFvMesh.H"
#include "cellSet.H"
#include "fvOptions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "createFvOptions.H"
// create cfdemCloud
#include "readGravitationalAcceleration.H"
cfdemCloudIB particleCloud(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
//=== dyM ===================
interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
mesh.update(); //dyM
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
// do particle stuff
Info << "- evolve()" << endl;
particleCloud.evolve(Us);
volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
// Pressure-velocity PISO corrector
{
MRF.correctBoundaryVelocity(U);
// Momentum predictor
#include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
#include "pEqn.H"
}
}
laminarTransport.correct();
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

143
applications/solvers/cfdemSolverIBContinuousForcing/createFields.H Executable file
View File

@ -0,0 +1,143 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading physical velocity field U" << endl;
Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading particle velocity field Us\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading the penalization factor field lambda\n" << endl;
volScalarField lambda
(
IOobject
(
"lambda",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// drag law modelling
//========================
Info<< "\nCreating dummy density field rho = 1\n" << endl;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
);
Info<< "Reading field phiIB\n" << endl;
volScalarField phiIB
(
IOobject
(
"phiIB",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//mod by alice
Info<< "Reading field voidfraction\n" << endl;
volScalarField voidfraction
(
IOobject
(
"voidfraction",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
//=== dyM ===================
Info<< "Reading field interFace\n" << endl;
volScalarField interFace
(
IOobject
(
"interFace",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
//dimensionedScalar("0", dimensionSet(0, -1, 0, 0, 0), 0.0)
dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 0.0)
);
//===========================
#include "createMRF.H"

35
applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H Normal file
View File

@ -0,0 +1,35 @@
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf(fvc::interpolate(rUA));
U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ rUAf*fvc::ddtCorr(U, phi); // Is there additional flux term due to the particle presence?
MRF.makeRelative(phi);
adjustPhi(phi, U, p);
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
#include "continuityErrs.H"
U -= rUA*fvc::grad(p); // should we add a pressure correction?
U.correctBoundaryConditions();
fvOptions.correct(U);

2
applications/solvers/cfdemSolverMultiphaseScalar/Allwclean
View File

@ -2,7 +2,7 @@
cd ${0%/*} || exit 1 # Run from this directory
set -x
wclean libso multiphaseMixture
wclean libso multiphaseMixtureScalar
wclean
#------------------------------------------------------------------------------

2
applications/solvers/cfdemSolverMultiphaseScalar/Allwmake
View File

@ -6,7 +6,7 @@ targetType=libso
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
set -x
wmake $targetType multiphaseMixture
wmake $targetType multiphaseMixtureScalar
wmake
#------------------------------------------------------------------------------

2
applications/solvers/cfdemSolverMultiphaseScalar/Make/options
View File

@ -6,7 +6,7 @@ include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
EXE_INC = \
$(PFLAGS) \
-I$(CFDEM_OFVERSION_DIR) \
-ImultiphaseMixture/lnInclude \
-ImultiphaseMixtureScalar/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

2
applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C
View File

@ -30,7 +30,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "multiphaseMixture.H"
#include "multiphaseMixtureScalar.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"

2
applications/solvers/cfdemSolverMultiphaseScalar/createFields.H
View File

@ -88,7 +88,7 @@ surfaceScalarField phi
linearInterpolate(U*voidfraction) & mesh.Sf()
);
multiphaseMixture mixture(U, phi, voidfraction);
multiphaseMixtureScalar mixture(U, phi, voidfraction);
// Need to store rho for ddt(rho, U)
volScalarField rho

2
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/Make/files
View File

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

0
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/Make/options
View File

0
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
View File

8
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/alphaContactAngle/alphaContactAngleFvPatchScalarField.H
View File

@ -26,7 +26,7 @@ Class
Description
Contact-angle boundary condition for multi-phase interface-capturing
simulations. Used in conjuction with multiphaseMixture.
simulations. Used in conjuction with multiphaseMixtureScalar.
SourceFiles
alphaContactAngleFvPatchScalarField.C
@ -37,7 +37,7 @@ SourceFiles
#define alphaContactAngleFvPatchScalarField_H
#include "zeroGradientFvPatchFields.H"
#include "multiphaseMixture.H"
#include "multiphaseMixtureScalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -117,8 +117,8 @@ public:
typedef HashTable
<
interfaceThetaProps,
multiphaseMixture::interfacePair,
multiphaseMixture::interfacePair::hash
multiphaseMixtureScalar::interfacePair,
multiphaseMixtureScalar::interfacePair::hash
> thetaPropsTable;

58
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.C
View File

@ -18,7 +18,7 @@ License
\*---------------------------------------------------------------------------*/
#include "multiphaseMixture.H"
#include "multiphaseMixtureScalar.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "Time.H"
#include "subCycle.H"
@ -31,13 +31,13 @@ License
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
const Foam::scalar Foam::multiphaseMixture::convertToRad =
const Foam::scalar Foam::multiphaseMixtureScalar::convertToRad =
Foam::constant::mathematical::pi/180.0;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::multiphaseMixture::calcAlphas()
void Foam::multiphaseMixtureScalar::calcAlphas()
{
scalar level = 0.0;
alphas_ == 0.0;
@ -51,7 +51,7 @@ void Foam::multiphaseMixture::calcAlphas()
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::calcNu() const
Foam::multiphaseMixtureScalar::calcNu() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -74,7 +74,7 @@ Foam::multiphaseMixture::calcNu() const
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::calcStf() const
Foam::multiphaseMixtureScalar::calcStf() const
{
tmp<surfaceScalarField> tstf
(
@ -134,7 +134,7 @@ Foam::multiphaseMixture::calcStf() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::multiphaseMixture::multiphaseMixture
Foam::multiphaseMixtureScalar::multiphaseMixtureScalar
(
const volVectorField& U,
const surfaceScalarField& phi,
@ -230,7 +230,7 @@ Foam::multiphaseMixture::multiphaseMixture
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::rho() const
Foam::multiphaseMixtureScalar::rho() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -247,7 +247,7 @@ Foam::multiphaseMixture::rho() const
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::rho(const label patchi) const
Foam::multiphaseMixtureScalar::rho(const label patchi) const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -264,9 +264,9 @@ Foam::multiphaseMixture::rho(const label patchi) const
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::mu() const
Foam::multiphaseMixtureScalar::mu() const
{
Info << "In multiphasemixture mu()" << endl;
Info << "In multiphaseMixtureScalar mu()" << endl;
return rho()*nu();
// PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -283,7 +283,7 @@ Foam::multiphaseMixture::mu() const
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::mu(const label patchi) const
Foam::multiphaseMixtureScalar::mu(const label patchi) const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -306,7 +306,7 @@ Foam::multiphaseMixture::mu(const label patchi) const
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::muf() const
Foam::multiphaseMixtureScalar::muf() const
{
return nuf()*fvc::interpolate(rho());
@ -327,13 +327,13 @@ Foam::multiphaseMixture::muf() const
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nu() const
Foam::multiphaseMixtureScalar::nu() const
{
return nu_;
}
Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::nu(const label patchi) const
Foam::multiphaseMixtureScalar::nu(const label patchi) const
{
//return nu_.boundaryField()[patchi];
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -355,7 +355,7 @@ Foam::multiphaseMixture::nu(const label patchi) const
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::nuf() const
Foam::multiphaseMixtureScalar::nuf() const
{
//return muf()/fvc::interpolate(rho());
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -374,7 +374,7 @@ Foam::multiphaseMixture::nuf() const
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cp() const
Foam::multiphaseMixtureScalar::Cp() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -395,7 +395,7 @@ Foam::multiphaseMixture::Cp() const
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::kf() const
Foam::multiphaseMixtureScalar::kf() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -417,7 +417,7 @@ Foam::multiphaseMixture::kf() const
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::D() const
Foam::multiphaseMixtureScalar::D() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -439,7 +439,7 @@ Foam::multiphaseMixture::D() const
}
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cs() const
Foam::multiphaseMixtureScalar::Cs() const
{
PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -456,7 +456,7 @@ Foam::multiphaseMixture::Cs() const
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::diffusionCorrection() const
Foam::multiphaseMixtureScalar::diffusionCorrection() const
{
surfaceScalarField numerator
@ -517,7 +517,7 @@ Foam::multiphaseMixture::diffusionCorrection() const
return correction;
}
void Foam::multiphaseMixture::solve()
void Foam::multiphaseMixtureScalar::solve()
{
correct();
@ -570,7 +570,7 @@ void Foam::multiphaseMixture::solve()
}
void Foam::multiphaseMixture::correct()
void Foam::multiphaseMixtureScalar::correct()
{
forAllIter(PtrDictionary<phase>, phases_, iter)
{
@ -579,7 +579,7 @@ void Foam::multiphaseMixture::correct()
}
Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureScalar::nHatfv
(
const volScalarField& alpha1,
const volScalarField& alpha2
@ -605,7 +605,7 @@ Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
}
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureScalar::nHatf
(
const volScalarField& alpha1,
const volScalarField& alpha2
@ -622,7 +622,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
// 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
void Foam::multiphaseMixtureScalar::correctContactAngle
(
const phase& alpha1,
const phase& alpha2,
@ -726,7 +726,7 @@ void Foam::multiphaseMixture::correctContactAngle
}
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureScalar::K
(
const phase& alpha1,
const phase& alpha2
@ -742,7 +742,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nearInterface() const
Foam::multiphaseMixtureScalar::nearInterface() const
{
tmp<volScalarField> tnearInt
(
@ -768,7 +768,7 @@ Foam::multiphaseMixture::nearInterface() const
}
void Foam::multiphaseMixture::solveAlphas
void Foam::multiphaseMixtureScalar::solveAlphas
(
const scalar cAlpha
)
@ -901,7 +901,7 @@ void Foam::multiphaseMixture::solveAlphas
}
bool Foam::multiphaseMixture::read()
bool Foam::multiphaseMixtureScalar::read()
{
if (transportModel::read())
{

18
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.H
View File

@ -17,10 +17,10 @@ License
Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
Class
multiphaseMixture
multiphaseMixtureScalar
Description
This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
This class is based on the OpenFOAM(R) Foam::multiphaseMixtureScalar 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
@ -33,11 +33,11 @@ Description
between each phase-pair.
SourceFiles
multiphaseMixture.C
multiphaseMixtureScalar.C
\*---------------------------------------------------------------------------*/
#ifndef multiphaseMixture_H
#define multiphaseMixture_H
#ifndef multiphaseMixtureScalar_H
#define multiphaseMixtureScalar_H
#include "incompressible/transportModel/transportModel.H"
#include "IOdictionary.H"
@ -52,10 +52,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class multiphaseMixture Declaration
Class multiphaseMixtureScalar Declaration
\*---------------------------------------------------------------------------*/
class multiphaseMixture
class multiphaseMixtureScalar
:
public IOdictionary,
public transportModel
@ -191,7 +191,7 @@ public:
// Constructors
//- Construct from components
multiphaseMixture
multiphaseMixtureScalar
(
const volVectorField& U,
const surfaceScalarField& phi,
@ -200,7 +200,7 @@ public:
//- Destructor
virtual ~multiphaseMixture()
virtual ~multiphaseMixtureScalar()
{}

0
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/phase/phase.C
View File

2
applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H → applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/phase/phase.H
View File

@ -26,7 +26,7 @@ Class
Description
Single incompressible phase derived from the phase-fraction.
Used as part of the multiPhaseMixture for interface-capturing multi-phase
Used as part of the multiphaseMixtureScalar for interface-capturing multi-phase
simulations.
SourceFiles

3
applications/solvers/cfdemSolverPisoFreeStreaming/Make/options
View File

@ -12,7 +12,8 @@ EXE_INC = \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-I$(LIB_SRC)/sampling/lnInclude
-I$(LIB_SRC)/sampling/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\

4
applications/solvers/cfdemSolverRhoPimple/EEqn.H
View File

@ -51,8 +51,4 @@
thermo.correct();
Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
particleCloud.clockM().start(31,"energySolve");
particleCloud.solve();
particleCloud.clockM().stop("energySolve");
}

4
applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
View File

@ -55,8 +55,4 @@ fvScalarMatrix EEqn
Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
Info << "he max/min : " << max(he).value() << " " << min(he).value() << endl;
particleCloud.clockM().start(31,"energySolve");
particleCloud.solve();
particleCloud.clockM().stop("energySolve");
}

4
applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
View File

@ -60,13 +60,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
if (propagateInertSpecie)
{
if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL);
if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
forAll(Y,i)
{
if (i!=inertIndex)
{
volScalarField& Yi = Y[i];
Yi = Yi/(Yt+VSMALL);
Yi = Yi/(Yt+ROOTVSMALL);
}
}
}

14
applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
View File

@ -303,18 +303,4 @@
mesh,
dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
);
volScalarField dSauter
(
IOobject
(
"dSauter",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
);
//===============================

11
applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H → applications/solvers/rcfdemSolverCoupledHeattransfer/TEqn.H
View File

@ -1,4 +1,4 @@
volScalarField rhoeps = rhoRec*voidfractionRec;
rhoeps = rhoRec*voidfractionRec;
particleCloud.energyContributions(Qsource);
@ -23,10 +23,17 @@
fvOptions(rhoeps, T) // no fvOptions support yet
);
fvOptions.constrain(TEqn); // no fvOptions support yet
// fvOptions.constrain(TEqn); // no fvOptions support yet
TEqn.relax();
TEqn.solve();
T = max(T, TMin);
T = min(T, TMax);
Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
particleCloud.clockM().start(31,"postFlow");
counter++;

25
applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H
View File

@ -73,6 +73,9 @@
dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
);
volScalarField rhoeps("rhoeps", rhoRec*voidfractionRec);
rhoeps.oldTime(); // switch on saving old time
// heat transfer fields
Info << "\nCreating heat transfer fields.\n" << endl;
@ -228,3 +231,25 @@
)
);
weightDict.add("weights",scalarList(1,1.0));
dimensionedScalar TMax
(
dimensionedScalar::lookupOrDefault
(
"TMax",
transportProps,
dimTemperature,
GREAT
)
);
dimensionedScalar TMin
(
dimensionedScalar::lookupOrDefault
(
"TMin",
transportProps,
dimTemperature,
0.0
)
);

2
applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C
View File

@ -98,7 +98,7 @@ int main(int argc, char *argv[])
particleCloud.clockM().start(26,"Flow");
#include "updateRho.H"
#include "TEqImp.H"
#include "TEqn.H"
particleCloud.clockM().stop("Flow");
stepCounter++;

3
applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H
View File

@ -1 +1,2 @@
rhoRec = pRec / (T * R);
dimensionedScalar Tave = T.weightedAverage(voidfractionRec);
rhoRec = pRec / (Tave * R);

3
applications/solvers/rcfdemSolverForcedTracers/Make/options
View File

@ -12,6 +12,7 @@ EXE_INC = \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\
@ -24,4 +25,4 @@ EXE_LIBS = \
-lfvOptions \
-l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS)
$(CFDEM_ADD_LIBS)

5
applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
View File

@ -60,9 +60,4 @@
thermo.correct();
Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
particleCloud.clockM().start(31,"energySolve");
particleCloud.solve();
particleCloud.clockM().stop("energySolve");
}

5
applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H
View File

@ -59,9 +59,4 @@
thermo.correct();
Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
particleCloud.clockM().start(31,"energySolve");
particleCloud.solve();
particleCloud.clockM().stop("energySolve");
}

1
applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options
View File

@ -27,6 +27,7 @@ EXE_INC = \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\

1
applications/solvers/rctfSpeciesTransport/Make/options
View File

@ -12,6 +12,7 @@ EXE_INC = \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-Wno-deprecated-copy
EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\

171
applications/utilities/displacementField/displacementField.C
View File

@ -37,11 +37,19 @@ Application
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void findPairs(labelList &, labelList &, labelPairList &);
void findPairsUnordered(labelList &, labelList &, labelPairList &);
void fillEmptyCells(fvMesh &, label, label, labelList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
void nearestNeighborCells(fvMesh &, label, label, label, labelList &, labelList &);
void normalizeFields(labelList &, volVectorField &, volVectorField &);
void readDump(std::string, labelList &, vectorList &);
void fillEmptyCells(fvMesh &, label, label, scalarList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
void nearestNeighborCells(fvMesh &, label, label, label, scalarList &, labelList &);
void normalizeFields(scalarList &, volVectorField &, volVectorField &);
void readDump(std::string, labelList &, scalarList &, vectorList &);
scalar weightFun(scalar);
label maxNumParticles = 1000000;
scalar minVol = 1e-12;
scalar Pi43 = 4.1888;
label posIndex = -1;
label posRad = -1;
label posX = -1;
label posY = -1;
label posZ = -1;
int main(int argc, char *argv[])
{
@ -88,6 +96,11 @@ int main(int argc, char *argv[])
label dumpIndexDisplacementIncrement(readLabel(displacementProperties.lookup("dumpIndexDisplacementIncrement")));
label nNeighMin(readLabel(displacementProperties.lookup("nNeighMin")));
label maxSearchLayers(displacementProperties.lookupOrDefault<label>("maxSearchLayers",0));
posIndex = readLabel(displacementProperties.lookup("posIndex"));
posRad = readLabel(displacementProperties.lookup("posRad"));
posX = readLabel(displacementProperties.lookup("posX"));
posY = readLabel(displacementProperties.lookup("posY"));
posZ = readLabel(displacementProperties.lookup("posZ"));
scalar timePerInputStep(readScalar(displacementProperties.lookup("timePerInputStep")));
scalar timePerDisplacementStep(readScalar(displacementProperties.lookup("timePerDisplacementStep")));
scalar startTime(readScalar(displacementProperties.lookup("startTime")));
@ -95,7 +108,7 @@ int main(int argc, char *argv[])
std::string fileext=string(displacementProperties.lookupOrDefault<string>("fileextension",""));
bool interpolate=bool(displacementProperties.lookupOrDefault<bool>("fillEmptyCells",true));
bool averageMode=bool(displacementProperties.lookupOrDefault<bool>("averageMode",false));
volVectorField defaultUs
(
IOobject
@ -110,11 +123,11 @@ int main(int argc, char *argv[])
dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
);
volVectorField defaultUsDirectedVariance
volVectorField defaultUsDirectedStdDev
(
IOobject
(
"defaultUDispDirectedVariance",
"defaultUDispDirectedStdDev",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
@ -172,11 +185,11 @@ int main(int argc, char *argv[])
dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
);
volVectorField UsDirectedVariance
volVectorField UsDirectedStdDev
(
IOobject
(
"UDispDirectedVariance",
"UDispDirectedStdDev",
runTime.timeName(),
mesh,
IOobject::NO_READ,
@ -186,7 +199,7 @@ int main(int argc, char *argv[])
dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
);
labelList particlesInCell(mesh.nCells(), 0);
scalarList particleVolInCell(mesh.nCells(), 0.0);
scalar currTime=startTime + thisProc * timePerInputStep;
label timeIndex=thisProc;
@ -196,6 +209,7 @@ int main(int argc, char *argv[])
runTime.setTime(currTime,timeIndex);
// read dump files and check which particle indices are present in both
labelList indices1, indices2;
scalarList radii1, radii2;
vectorList positions1, positions2;
std::stringstream ss;
@ -209,13 +223,13 @@ int main(int argc, char *argv[])
{
if (averageMode)
{
normalizeFields(particlesInCell, Us, UsDirectedVariance);
fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
Us /= timePerDisplacementStep;
UsDirectedVariance /= timePerDisplacementStep;
UsDirectedStdDev /= timePerDisplacementStep;
Us.write();
UsDirectedVariance.write();
UsDirectedStdDev.write();
}
break;
}
@ -225,8 +239,8 @@ int main(int argc, char *argv[])
Info << "\t" << filename2 << endl;
Info << "corresponding to time = " << currTime << "." << endl;
readDump(filename1, indices1, positions1);
readDump(filename2, indices2, positions2);
readDump(filename1, indices1, radii1, positions1);
readDump(filename2, indices2, radii2, positions2);
labelPairList pairs;
findPairs(indices1,indices2,pairs);
@ -234,15 +248,16 @@ int main(int argc, char *argv[])
// average particle displacements and their variance
Info << "Binning particle displacements on mesh." << endl;
vector position, displacement;
scalar radius, volume;
label line1, line2;
label cellI;
if (!averageMode)
{
Us *= 0.0;
UsDirectedVariance *= 0.0;
particlesInCell.clear();
particlesInCell.setSize(mesh.nCells(), 0);
UsDirectedStdDev *= 0.0;
particleVolInCell.clear();
particleVolInCell.setSize(mesh.nCells(), 0);
}
for (label partI = 0; partI < pairs.size(); partI++)
@ -250,27 +265,29 @@ int main(int argc, char *argv[])
line1 = pairs[partI].first();
line2 = pairs[partI].second();
position = positions1[line1];
displacement = positions2[line2] - positions1[line1];
cellI = mesh.findCell(position);
if (cellI < 0) continue;
particlesInCell[cellI] += 1;
Us[cellI] += displacement;
displacement = positions2[line2] - positions1[line1];
radius = radii1[line1];
volume = Pi43 * radius * radius * radius;
particleVolInCell[cellI] += volume;
Us[cellI] += displacement*volume;
for (label comp=0;comp<3;comp++)
{
UsDirectedVariance[cellI].component(comp) += displacement.component(comp)*displacement.component(comp);
UsDirectedStdDev[cellI].component(comp) += displacement.component(comp)*displacement.component(comp)*volume;
}
}
if (!averageMode)
{
normalizeFields(particlesInCell, Us, UsDirectedVariance);
fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
Us /= timePerDisplacementStep;
UsDirectedVariance /= timePerDisplacementStep;
UsDirectedStdDev /= timePerDisplacementStep;
Us.write();
UsDirectedVariance.write();
UsDirectedStdDev.write();
}
if (averageMode && monitorProbes)
@ -280,7 +297,7 @@ int main(int argc, char *argv[])
{
vector pos = probePoints[p];
label cellP = mesh.findCell(pos);
monitoringDataFile << " " << particlesInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedVariance[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
monitoringDataFile << " " << particleVolInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedStdDev[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
}
monitoringDataFile << endl;
}
@ -293,30 +310,69 @@ int main(int argc, char *argv[])
return 0;
}
void readDump(std::string filename, labelList &indices, vectorList &positions)
void readDump(std::string filename, labelList &indices, scalarList &radii, vectorList &positions)
{
#include <fstream>
const label leadingLines = 9;
label lineCounter = 0;
label partIndex;
scalar x, y, z;
scalar r = 1.0, x = 0.0, y = 0.0, z = 0.0;
indices.clear();
radii.clear();
positions.clear();
indices.setSize(maxNumParticles);
radii.setSize(maxNumParticles);
positions.setSize(maxNumParticles);
std::ifstream file(filename);
std::string str;
std::string word;
label wordcounter;
while (std::getline(file, str))
{
if (lineCounter >= leadingLines)
{
sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
indices.append(partIndex);
positions.append(vector(x,y,z));
std::istringstream ss(str);
wordcounter = 0;
while (ss >> word)
{
if (wordcounter == posIndex)
{
partIndex = stoi(word);
}
else if (wordcounter == posRad)
{
r = stod(word);
}
else if (wordcounter == posX)
{
x = stod(word);
}
else if (wordcounter == posY)
{
y = stod(word);
}
else if (wordcounter == posZ)
{
z = stod(word);
}
wordcounter++;
}
// sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
indices[lineCounter-leadingLines] = partIndex;
radii[lineCounter-leadingLines] = r;
positions[lineCounter-leadingLines] = vector(x,y,z);
}
lineCounter++;
}
label readLines = lineCounter - leadingLines;
indices.resize(readLines);
radii.resize(readLines);
positions.resize(readLines);
}
void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
@ -324,6 +380,10 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
// remove all entries from first list if they are not present in second list
// this assumes ordered entries
pairs.clear();
pairs.setSize(maxNumParticles);
label pairCounter = 0;
if (indices2.size() == 0) return;
for (label i=0;i<indices1.size();i++)
@ -339,18 +399,23 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
else if (indices2[jmid] < index1) j1 = jmid;
else
{
pairs.append(labelPair(i,jmid));
pairs[pairCounter]=labelPair(i,jmid);
pairCounter++;
break;
}
if (j2-j1 == 1) break;
}
}
pairs.resize(pairCounter);
Info << "findPairs: " << pairs.size() << " pairs found." << endl;
}
void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList &pairs)
{
// remove all entries from first list if they are not present in second list
pairs.clear();
pairs.setSize(maxNumParticles);
label pairCounter = 0;
for (label i=0;i<indices1.size();i++)
{
@ -358,15 +423,17 @@ void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList
{
if (indices1[i] == indices2[j])
{
pairs.append(labelPair(i,j));
pairs[pairCounter]=labelPair(i,j);
pairCounter++;
break;
}
}
}
pairs.resize(pairCounter);
Info << "findPairs: " << pairs.size() << " pairs found." << endl;
}
void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, volVectorField &Us, volVectorField& UsDirectedVariance,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedVariance, bool interpolate, scalar dt)
void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, volVectorField &Us, volVectorField& UsDirectedStdDev,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedStdDev, bool interpolate, scalar dt)
{
labelList neighborsWithValues;
scalar neighborSqrDistance;
@ -377,7 +444,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
Info << "Filling empty cells." << endl;
forAll(mesh.C(), cellI)
{
if (particlesInCell[cellI] > 0) continue;
if (particleVolInCell[cellI] > minVol) continue;
vector position = mesh.C()[cellI];
label outsideBox = 0;
@ -388,11 +455,11 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
if (outsideBox > 0 || !interpolate)
{
Us[cellI] = defaultUs[cellI]*dt;
UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
continue;
}
nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particlesInCell, neighborsWithValues);
nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particleVolInCell, neighborsWithValues);
weightSum = 0.0;
weights.clear();
for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
@ -406,13 +473,13 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
{
weight = weights[neighI]/weightSum;
Us[cellI] += weight*Us[neighborsWithValues[neighI]];
UsDirectedVariance[cellI] += weight*UsDirectedVariance[neighborsWithValues[neighI]];
UsDirectedStdDev[cellI] += weight*UsDirectedStdDev[neighborsWithValues[neighI]];
}
if (neighborsWithValues.size() == 0)
{
Us[cellI] = defaultUs[cellI]*dt;
UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
}
// make sure no particles are placed outside of domain
@ -429,7 +496,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
}
}
void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, labelList &neighborsWithValues)
void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, labelList &neighborsWithValues)
{
label numSearchLayers = 0;
std::set<label> neighbors;
@ -455,11 +522,11 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
{
newNeighbors.insert(adj);
neighbors.insert(adj);
if (particlesInCell[adj] > 0) neighborsWithValues.append(adj);
if (particleVolInCell[adj] > minVol) neighborsWithValues.append(adj);
}
}
}
numSearchLayers++;
if (numSearchLayers > maxSearchLayers && maxSearchLayers > 0) return;
@ -470,18 +537,18 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
}
}
void normalizeFields(labelList& particlesInCell, volVectorField& Us, volVectorField & UsDirectedVariance)
void normalizeFields(scalarList& particleVolInCell, volVectorField& Us, volVectorField & UsDirectedStdDev)
{
for (label cellJ = 0; cellJ<particlesInCell.size(); cellJ++)
for (label cellJ = 0; cellJ<particleVolInCell.size(); cellJ++)
{
if (particlesInCell[cellJ] > 0)
if (particleVolInCell[cellJ] > minVol)
{
Us[cellJ] /= particlesInCell[cellJ];
UsDirectedVariance[cellJ] /= particlesInCell[cellJ];
Us[cellJ] /= particleVolInCell[cellJ];
UsDirectedStdDev[cellJ] /= particleVolInCell[cellJ];
for (label comp=0;comp<3;comp++)
{
UsDirectedVariance[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
if (UsDirectedVariance[cellJ].component(comp) > 0) UsDirectedVariance[cellJ].component(comp) = Foam::sqrt(UsDirectedVariance[cellJ].component(comp));
UsDirectedStdDev[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
if (UsDirectedStdDev[cellJ].component(comp) > 0) UsDirectedStdDev[cellJ].component(comp) = Foam::sqrt(UsDirectedStdDev[cellJ].component(comp));
}
}
}

1
doc/CFDEMcoupling_models.txt
View File

@ -124,6 +124,7 @@ particleDeformation,
potentialRelaxation,
"surfaceTensionForce"_forceModel_surfaceTensionForce.html,
terminalVelocity,
"transferFluidProperties"_forceModel_transferFluidProperties.html,
turbulentDispersion,
turbulentVelocityFluctuations,
"virtualMassForce"_forceModel_virtualMassForce.html,

1
doc/CFDEMcoupling_solvers.txt
View File

@ -10,6 +10,7 @@
This section lists all CFDEMcoupling solvers alphabetically.
"cfdemSolverIB"_cfdemSolverIB.html,
"cfdemSolverIBContinuousForcing"_cfdemSolverIBContinuousForcing.html,
"cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
"cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html,
"cfdemSolverPiso"_cfdemSolverPiso.html,

8
doc/IOModel_basicIO.txt
View File

@ -17,7 +17,8 @@ IOModel basicIO; :pre
[Examples:]
IOModel basicIO;
serialOutput; :pre
serialOutput;
cartesianOutput; :pre
[Description:]
@ -29,6 +30,11 @@ Using the keyword {serialOutput;} in the
the IO is serial to the directory {$casePath/CFD/lagrangian}. In this case
only the data on processor 0 is written!
Using the keyword {cartesianOutput;} in the
"couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties dictionary,
particle positions will be written as cartesian coordinates instead of
barycentric coordinates.
Data is written every write time of the CFD simulation.
[Restrictions:]

82
doc/cfdemSolverIBContinuousForcing.txt Normal file
View File

@ -0,0 +1,82 @@
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
cfdemSolverIBContinuousForcing command :h3
[Description:]
"cfdemSolverIBContinuousForcing" is a coupled CFD-DEM solver using CFDEMcoupling,
an open-source parallel coupled CFD-DEM framework, for calculating the dynamics
between immersed bodies and the surrounding fluid. Being an implementation of a
continuous forcing immersed boundary method it allows tackling problems where
the body diameter exceeds the maximal size of a fluid cell.
<!-- HTML_ONLY -->
Using the toolbox of OpenFOAM&reg;(*) the governing equations of the fluid are
computed and the corrections of velocity and pressure field with respect to the
body-movement information, gained from LIGGGHTS, are incorporated.
<!-- END_HTML_ONLY -->
<!-- RST
Using the toolbox of OpenFOAM\ |reg|\ (*) the governing equations of the fluid are
computed and the corrections of velocity and pressure field with respect to the
body-movement information, gained from LIGGGHTS, are incorporated.
.. |reg| unicode:: U+000AE .. REGISTERED SIGN
END_RST -->
The code of this solver was contributed by A.N. Balachandran Nair, JKU. For more
details, see "Balachandran Nair et al. (2021)"_#BalachandranNair2021
[Use:]
The solver is realized within the open-source framework CFDEMcoupling. Just as
for the unresolved CFD-DEM solver cfdemSolverPiso the file
CFD/constant/couplingProperties contains information about the settings for the
different models. While IOmodel, DataExchangeModel etc. are applicable for all
CFDEMcoupling-solvers, special locate-, force- and void fraction models were
designed for this solver:
"engineSearchIB"_locateModel_engineSearchIB.html,
"ArchimedesIB"_forceModel_ArchimedesIB.html,
"ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
"IBVoidfraction"_voidFractionModel_IBVoidFraction.html
:line
:link(BalachandranNair2021)
[(Balachandran Nair, 2021)] Balachandran Nair, A.N., Pirker, S. and Saeedipour, M.,
"Resolved CFD-DEM simulation of blood flow with a reduced-order RBC model",
Comp. Part. Mech. (2021)
:line
<!-- HTML_ONLY -->
NOTE:
(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com, and owner of the
OPENFOAM&reg; and OpenCFD&reg; trade marks.
OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com.
<!-- END_HTML_ONLY -->
<!-- RST
.. note::
(*) This offering is not approved or endorsed by OpenCFD Limited, producer
and distributor of the OpenFOAM software via www.openfoam.com, and owner of
the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
distributor of the OpenFOAM software via www.openfoam.com.
.. |reg| unicode:: U+000AE .. REGISTERED SIGN
END_RST -->

10
doc/forceModel_ParmarBassetForce.txt
View File

@ -23,17 +23,19 @@ ParmarBassetForceProps
nIntegral scalar1;
discretisationOrder scalar2;
treatForceExplicit switch1;
interpolation switch2;
verbose switch2;
interpolation switch3;
smoothingModel "smoothingModel";
\} :pre
{U} = name of the finite volume fluid velocity field :ulb,l
{Us} = name of the finite volume cell averaged particle velocity field :l
{scalar1} = number of timesteps considered in the near history :l
{scalar2} = (1 or 2) discretisation order of the far history differential equations :l
{scalar1} = number of timesteps considered in the near history (integer > 1):l
{scalar2} = (optional, default 1) discretisation order of the far history differential equations (1 or 2) :l
{switch1} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{smoothingModel} = name of smoothing model for the dU/dt field :l
{switch3} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
{smoothingModel} = name of smoothing model for the dUrel/dt field :l
:ule
[Examples:]

2
doc/forceModel_ShirgaonkarIB.txt
View File

@ -21,12 +21,14 @@ ShirgaonkarIBProps
velFieldName "U";
pressureFieldName "pressure";
twoDimensional;
useTorque;
treatForceExplicit switch1;
\} :pre
{U} = name of the finite volume fluid velocity field :ulb,l
{pressure} = name of the finite volume pressure field :l
{twoDimensional} = optional keyword for conducting a two dimensional calculation :l
{useTorque} = optional keyword for calculating particle torque :l
{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
:ule

42
doc/forceModel_transferFluidProperties.txt Normal file
View File

@ -0,0 +1,42 @@
"CFDEMproject Website"_lws - "Main Page"_main :c
:link(lws,http://www.cfdem.com)
:link(main,CFDEMcoupling_Manual.html)
:line
forceModel transferFluidProperties command :h3
[Syntax:]
Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
dictionary.
forceModels
(
transferFluidProperties
);
transferFluidPropertiesProps
\{
verbose switch1;
interpolation switch2;
\} :pre
{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :ulb,l
{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
:ule
[Description:]
This "force model" does not influence the particles or the flow - it transfer to fluid density and (dynamic)
viscosity from OpenFOAM to LIGGGHTS.
[Restrictions:]
This model requires {fix cfd/coupling/fluidproperties} to work.
[Related commands:]
"forceModel"_forceModel.html

20
doc/smoothingModel_constDiffSmoothing.txt
View File

@ -15,17 +15,22 @@ dictionary.
smoothingModel constDiffSmoothing;
constDiffSmoothingProps
\{
lowerLimit number1;
upperLimit number2;
smoothingLength lengthScale;
smoothingLengthReferenceField lengthScaleRefField;
lowerLimit number1;
upperLimit number2;
smoothingLength lengthScale;
smoothingLengthReference lengthScaleRef;
smoothingLengthFieldName fieldName1;
smoothingLengthReferenceFieldName fieldName2;
verbose;
\} :pre
{number1} = scalar fields will be bound to this lower value :ulb,l
{number2} = scalar fields will be bound to this upper value :l
{lengthScale} = length scale over which the exchange fields will be smoothed out :l
{lengthScaleRefField} = length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
{lengthScaleRef} = (optional) length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
{fieldName1} = (optional) name of scalar field to be used as local smoothing length. :l
{fieldName2} = (optional) name of scalar field to be used as local smoothing length for reference fields. :l
{verbose} = (optional, default false) flag for debugging output :l
:ule
@ -51,6 +56,11 @@ which these reference fields are not specified. Values calculated in the cells
e.g. the average particle velocity, which are not specified in all cells in case
the flow is rather dilute.
Alternative to {smoothingLength} and {smoothingLengthReference},
{smoothingLengthFieldName} and/or {smoothingLengthReferenceFieldName} can be used
to define spatial variation of the smoothing lengths. Either the scalar or field
options must be used, giving both will result in errors.
[Restrictions:]
This model is tested in a limited number of flow situations.

2
etc/bashrc
View File

@ -17,7 +17,7 @@
#------------------------------------------------------------------------------
export CFDEM_PROJECT=CFDEM
export CFDEM_VERSION=21.11
export CFDEM_VERSION=23.02
################################################################################
# USER EDITABLE PART: Changes made here may be lost with the next upgrade

2
etc/cshrc
View File

@ -15,7 +15,7 @@
#------------------------------------------------------------------------------
setenv CFDEM_PROJECT CFDEM
setenv CFDEM_VERSION 21.11
setenv CFDEM_VERSION 23.02
################################################################################
# USER EDITABLE PART: Changes made here may be lost with the next upgrade

2
etc/functions.sh
View File

@ -359,7 +359,7 @@ cleanCFDEMcase()
cd $casePath/CFD
cleanCase
rm -r $casePath/DEM/post/*
echo "dummyfile" >> $casePath/DEM/post/dummy
touch $casePath/DEM/post/.gitignore
cd $casePath
echo "done"
}

2
etc/library-list.txt
View File

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

1
etc/solver-list.txt
View File

@ -10,6 +10,7 @@ cfdemSolverPiso/dir
cfdemSolverRhoPimple/dir
cfdemSolverIB/dir
cfdemSolverPisoScalar/dir
cfdemSolverIBContinuousForcing/dir
cfdemSolverRhoPimpleChem/dir
cfdemSolverMultiphase/dir
cfdemSolverMultiphaseScalar/dir

2
etc/tutorial-list.txt
View File

@ -19,3 +19,5 @@ cfdemSolverIB/twoSpheresGlowinskiMPI/dir
cfdemSolverPisoScalar/packedBedTemp/dir
cfdemSolverPiso/ErgunTestCG/dir
cfdemSolverIB/falling_sphere_two_way_coupling/dir

1
src/lagrangian/cfdemParticle/Make/files
View File

@ -96,6 +96,7 @@ $(forceModels)/potentialRelaxation/potentialRelaxation.C
$(forceModels)/BeetstraDrag/BeetstraDrag.C
$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
$(forceModels)/dSauter/dSauter.C
$(forceModels)/transferFluidProperties/transferFluidProperties.C
$(forceModels)/Fines/Fines.C
$(forceModels)/Fines/FinesFields.C
$(forceModels)/Fines/FanningDynFines.C

4
src/lagrangian/cfdemParticle/cfdTools/versionInfo.H
View File

@ -34,8 +34,8 @@ Description
#ifndef versionInfo_H
#define versionInfo_H
word CFDEMversion="PFM 21.11";
word compatibleLIGGGHTSversion="PFM 21.11";
word CFDEMversion="PFM 23.02";
word compatibleLIGGGHTSversion="PFM 23.02";
word OFversion="6";
Info << "\nCFDEMcoupling version: " << CFDEMversion << endl;

14
src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C
View File

@ -274,6 +274,8 @@ bool cfdemCloudEnergy::evolve
void cfdemCloudEnergy::postFlow()
{
if (ignore()) return;
cfdemCloud::postFlow();
forAll(energyModel_, modeli)
{
@ -285,18 +287,6 @@ void cfdemCloudEnergy::postFlow()
}
}
void cfdemCloudEnergy::solve()
{
forAll(energyModel_, modeli)
{
energyModel_[modeli].solve();
}
forAll(massTransferModel_, modeli)
{
massTransferModel_[modeli].solve();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

4
src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.H
View File

@ -42,7 +42,7 @@ namespace Foam
class energyModel;
class massTransferModel;
class thermCondModel;
class diffCoeffModel;
class diffCoeffModel;
class chemistryModel;
/*---------------------------------------------------------------------------*\
Class cfdemCloudEnergy Declaration
@ -142,8 +142,6 @@ public:
void postFlow();
void solve();
};

59
src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C
View File

@ -36,7 +36,6 @@ Description
#include "locateModel.H"
#include "dataExchangeModel.H"
#include "IOModel.H"
#include <mpi.h>
#include "IOmanip.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -58,6 +57,7 @@ cfdemCloudIB::cfdemCloudIB
angularVelocities_(NULL),
pRefCell_(readLabel(mesh.solutionDict().subDict("PISO").lookup("pRefCell"))),
pRefValue_(readScalar(mesh.solutionDict().subDict("PISO").lookup("pRefValue"))),
DEMTorques_(NULL),
haveEvolvedOnce_(false),
skipLagrangeToEulerMapping_(false)
{
@ -75,6 +75,7 @@ cfdemCloudIB::cfdemCloudIB
cfdemCloudIB::~cfdemCloudIB()
{
dataExchangeM().destroy(angularVelocities_,3);
dataExchangeM().destroy(DEMTorques_,3);
}
@ -82,8 +83,7 @@ cfdemCloudIB::~cfdemCloudIB()
void cfdemCloudIB::getDEMdata()
{
cfdemCloud::getDEMdata();
Info << "=== cfdemCloudIB::getDEMdata() === particle rotation not considered in CFD" << endl;
//dataExchangeM().getData("omega","vector-atom",angularVelocities_);
dataExchangeM().getData("omega","vector-atom",angularVelocities_);
}
bool cfdemCloudIB::reAllocArrays()
@ -92,12 +92,24 @@ bool cfdemCloudIB::reAllocArrays()
{
// get arrays of new length
dataExchangeM().allocateArray(angularVelocities_,0.,3);
dataExchangeM().allocateArray(DEMTorques_,0.,3);
return true;
}
return false;
}
bool cfdemCloudIB::evolve()
void cfdemCloudIB::giveDEMdata()
{
cfdemCloud::giveDEMdata();
dataExchangeM().giveData("hdtorque","vector-atom",DEMTorques_);
}
inline double ** cfdemCloudIB::DEMTorques() const
{
return DEMTorques_;
}
bool cfdemCloudIB::evolve(volVectorField& Us)
{
numberOfParticlesChanged_ = false;
arraysReallocated_=false;
@ -140,6 +152,10 @@ bool cfdemCloudIB::evolve()
for (int i=0;i<nrForceModels();i++) forceM(i).setForce();
if(verbose_) Info << "setForce done." << endl;
// set particle velocity field
if(verbose_) Info << "- setVelocity(velocities_)" << endl;
calcForcingTerm(Us); // only needed for continuous forcing
// write DEM data
if(verbose_) Info << " -giveDEMdata()" << endl;
giveDEMdata();
@ -148,7 +164,8 @@ bool cfdemCloudIB::evolve()
haveEvolvedOnce_=true;
}
Info << "evolve done." << endl;
if(verbose_) Info << "evolve done." << endl;
//if(verbose_) #include "debugInfo.H";
@ -158,6 +175,34 @@ bool cfdemCloudIB::evolve()
return doCouple;
}
void cfdemCloudIB::calcForcingTerm(volVectorField& Us)
{
label cell = 0;
vector uP(0,0,0);
vector rVec(0,0,0);
vector velRot(0,0,0);
vector angVel(0,0,0);
Us.primitiveFieldRef() = vector::zero;
for (int index = 0; index < numberOfParticles(); ++index)
{
for (int subCell = 0; subCell < voidFractionM().cellsPerParticle()[index][0]; subCell++)
{
cell = cellIDs()[index][subCell];
if (cell >=0)
{
// calc particle velocity
for (int i=0;i<3;i++) rVec[i]=Us.mesh().C()[cell][i]-position(index)[i];
for (int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i];
velRot=angVel^rVec;
for (int i=0;i<3;i++) uP[i] = velocities()[index][i]+velRot[i];
Us[cell] = (1-voidfractions_[index][subCell])*uP;
}
}
}
Us.correctBoundaryConditions();
}
void cfdemCloudIB::calcVelocityCorrection
(
volScalarField& p,
@ -190,7 +235,7 @@ void cfdemCloudIB::calcVelocityCorrection
for(int i=0;i<3;i++) uParticle[i] = velocities()[index][i]+velRot[i];
// impose field velocity
U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
U[cellI]= (1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
}
}
//}
@ -213,7 +258,7 @@ void cfdemCloudIB::calcVelocityCorrection
U.correctBoundaryConditions();
// correct the pressure as well
p=p+phiIB/U.mesh().time().deltaT(); // do we have to account for rho here?
p=p+phiIB/U.mesh().time().deltaT(); // no need to account for rho since p=(p/rho) in OF
p.correctBoundaryConditions();
if (couplingProperties_.found("checkinterface"))

11
src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H
View File

@ -62,6 +62,8 @@ protected:
label pRefCell_;
scalar pRefValue_;
double **DEMTorques_;
bool haveEvolvedOnce_;
bool skipLagrangeToEulerMapping_;
@ -85,7 +87,13 @@ public:
bool reAllocArrays();
bool evolve();
void giveDEMdata();
inline double ** DEMTorques() const;
bool evolve(volVectorField&);
void calcForcingTerm(volVectorField&);
void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model
@ -96,7 +104,6 @@ public:
{
return angularVelocities_;
}
};

43
src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C
View File

@ -125,6 +125,14 @@ void IOModel::streamDataToPath(const fileName& path, const double* const* array,
else fileStream << "{}" << "\n";
return;
}
#if OPENFOAM_VERSION_MAJOR > 4
else if (barycentricOutput_ && type == "position")
{
// Force construction of face-diagonal decomposition before construction
// of particle which uses parallel transfers.
(void)particleCloud_.mesh().tetBasePtIs();
}
#endif
fileStream << token::BEGIN_LIST << nl;
@ -139,19 +147,24 @@ void IOModel::streamDataToPath(const fileName& path, const double* const* array,
}
else if (type == "position")
{
#if OPENFOAM_VERSION_MAJOR < 5
fileStream << "( " << array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << nl;
#else
particle part
(
particleCloud_.mesh(),
vector(array[index][0],array[index][1],array[index][2]),
cellIDs[index][0]
);
#if OPENFOAM_VERSION_MAJOR > 4
if (barycentricOutput_)
{
particle part
(
particleCloud_.mesh(),
vector(array[index][0],array[index][1],array[index][2]),
cellIDs[index][0]
);
part.writePosition(fileStream);
fileStream << nl;
part.writePosition(fileStream);
fileStream << nl;
}
else
#endif
{
fileStream << "( " << array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << nl;
}
}
else if (type == "vector")
{
@ -179,7 +192,8 @@ IOModel::IOModel
dict_(dict),
particleCloud_(sm),
time_(sm.mesh().time()),
parOutput_(true)
parOutput_(true),
barycentricOutput_(true)
{
if (
particleCloud_.dataExchangeM().type()=="oneWayVTK" ||
@ -189,6 +203,11 @@ IOModel::IOModel
parOutput_ = false;
Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
}
if (dict_.found("cartesianOutput"))
{
barycentricOutput_ = false;
}
}

2
src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.H
View File

@ -65,6 +65,8 @@ protected:
bool parOutput_;
bool barycentricOutput_;
public:
//- Runtime type information

2
src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
View File

@ -246,8 +246,6 @@ void diffusionCoefficient::execute()
double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_);
}
Info << "give data done" << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

31
src/lagrangian/cfdemParticle/subModels/chemistryModel/initMultiLayers/initMultiLayers.C
View File

@ -57,8 +57,10 @@ initMultiLayers::initMultiLayers
mesh_(sm.mesh()),
verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
maxNumLayers_(0),
maxNumParticlesPerType_(propsDict_.lookupOrDefault<label>("maxNumParticlesPerType",1000000)),
numLayers_(propsDict_.lookupOrDefault<labelList>("numLayers",labelList(1,-1))),
partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
searchLayers_(propsDict_.lookupOrDefault<labelList>("searchLayers",labelList(1,-1))),
relRadiiRegName_(typeName + "relRadii"),
filepath_(string(propsDict_.lookup("filepath"))),
initialized_(false)
@ -107,7 +109,7 @@ bool initMultiLayers::init()
}
Info << "\nReading" << endl;
Info << "\t" << filename << endl;
readDump(filename, type, indices, positions, relradii);
label numLines = readDump(filename, type, indices, positions, relradii);
Info << "Binning particle displacements on mesh for type " << type << endl;
@ -144,7 +146,7 @@ bool initMultiLayers::init()
);
label cellI = -1;
for (label lineI = 0; lineI < indices.size(); lineI++)
for (label lineI = 0; lineI < numLines; lineI++)
{
position = positions[lineI];
relrad = relradii[lineI];
@ -178,8 +180,8 @@ bool initMultiLayers::init()
if(cellK >= 0 && partType == type)
{
// look for the nearest occupied cell
cellKoccupied = findNearestCellWithValue(cellK, particlesInCell);
listIndex = getListIndex(partType);
cellKoccupied = findNearestCellWithValue(cellK, particlesInCell,searchLayers_[listIndex]);
if (cellKoccupied >= 0)
{
relRadiiLoopVar = relRadiiField[cellKoccupied];
@ -215,7 +217,7 @@ void initMultiLayers::execute()
}
}
label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &particlesInCell) const
label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &particlesInCell, label searchLayers) const
{
if (particlesInCell[refCell] > 0.5) return refCell;
@ -226,6 +228,8 @@ label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &p
neighbors.insert(refCell);
recentNeighbors.insert(refCell);
label layersSearched = 0;
while(true)
{
for (std::set<label>::iterator it=recentNeighbors.begin(); it!=recentNeighbors.end(); ++it)
@ -250,6 +254,8 @@ label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &p
recentNeighbors.clear();
recentNeighbors = newNeighbors;
newNeighbors.clear();
layersSearched++;
if (layersSearched >= searchLayers && searchLayers > 0) return -1;
}
}
@ -264,7 +270,7 @@ label initMultiLayers::getListIndex(label testElement) const
return -1;
}
void initMultiLayers::readDump(std::string filename, label type, labelList &indices, vectorList &positions, vectorList &relradii)
label initMultiLayers::readDump(std::string filename, label type, labelList &indices, vectorList &positions, vectorList &relradii)
{
#include <fstream>
@ -281,6 +287,10 @@ void initMultiLayers::readDump(std::string filename, label type, labelList &indi
positions.clear();
relradii.clear();
indices.setSize(maxNumParticlesPerType_);
positions.setSize(maxNumParticlesPerType_);
relradii.setSize(maxNumParticlesPerType_);
std::ifstream file(filename);
std::string str;
while (std::getline(file, str))
@ -292,12 +302,17 @@ void initMultiLayers::readDump(std::string filename, label type, labelList &indi
{
FatalError<< "Particle of type " << partType << " detected in " << filename << abort(FatalError);
}
indices.append(partIndex);
positions.append(vector(x,y,z));
relradii.append(vector(r1,r2,r3));
indices[lineCounter-leadingLines] = partIndex;
positions[lineCounter-leadingLines] = vector(x,y,z);
relradii[lineCounter-leadingLines] = vector(r1,r2,r3);
}
lineCounter++;
if (lineCounter == maxNumParticlesPerType_) break;
}
label readLines = lineCounter - leadingLines;
return readLines;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

8
src/lagrangian/cfdemParticle/subModels/chemistryModel/initMultiLayers/initMultiLayers.H
View File

@ -59,11 +59,15 @@ private:
bool verbose_;
label maxNumLayers_;
label maxNumParticlesPerType_;
labelList numLayers_;
labelList partTypes_;
labelList searchLayers_;
vector defaultRelRadii_;
const word relRadiiRegName_;
@ -74,9 +78,9 @@ private:
bool init();
void readDump(std::string, label, labelList &, vectorList &, vectorList &);
label readDump(std::string, label, labelList &, vectorList &, vectorList &);
label findNearestCellWithValue(label, volScalarField &) const;
label findNearestCellWithValue(label, volScalarField &, label) const;
label getListIndex(label) const;

2
src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
View File

@ -174,8 +174,6 @@ void massTransferCoeff::execute()
// give DEM data
particleCloud_.dataExchangeM().giveData(partNuName_, "scalar-atom", partNu_);
particleCloud_.dataExchangeM().giveData(partReName_, "scalar-atom", partRe_);
Info << "give data done" << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

2
src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
View File

@ -169,8 +169,6 @@ void reactantPerParticle::execute()
// give DEM data
particleCloud_.dataExchangeM().giveData(partReactantName_, "scalar-atom", reactantPerParticle_);
Info << "give data done" << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

127
src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.C
View File

@ -100,6 +100,7 @@ void One2One::setup
comm_
);
// count number of receives
ncollected_ = 0;
int nrequests = 0;
for (int i = 0; i < nsrc_procs_; i++)
@ -115,6 +116,18 @@ void One2One::setup
}
}
// count number of sends
if (nlocal_ > 0)
{
for (int i = 0; i < ndst_procs_; i++)
{
if (dst_procs_[i] != me_)
{
nrequests++;
}
}
}
if (nrequests > 0)
{
request_ = new MPI_Request[nrequests];
@ -143,23 +156,27 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
{
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
int tag = (src_procs_[i] << 16) | me_;
#ifdef O2ODEBUG
std::cout<< "[" << me_ << "]"
<< " RCV " << (i+1)
<< " of " << nsrc_procs_
<< " from: " << src_procs_[i]
<< " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]]
<< " data_length " << data_length
<< " offset " << offset
<< " tag " << tag
<< std::endl;
#endif
MPI_Irecv
(
&dst[offset],
natoms_[src_procs_[i]]*data_length,
wrap.mpi_type,
src_procs_[i],
MPI_ANY_TAG,
tag,
comm_,
&request_[requesti]
);
@ -168,16 +185,34 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
else // data is available on-proc
{
offset_local = offset;
#ifdef O2ODEBUG
std::cout<< "[" << me_ << "]"
<< " RCV " << (i+1)
<< " of " << nsrc_procs_
<< " from: " << src_procs_[i]
<< " skipped (self)"
<< " offset_local " << offset_local
<< std::endl;
#endif
}
}
#ifdef O2ODEBUG
else
{
std::cout<< "[" << me_ << "]"
<< " RCV " << (i+1)
<< " of " << nsrc_procs_
<< " from: " << src_procs_[i]
<< " skipped (no data)"
<< std::endl;
}
#endif
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?
// nonblocking sends
// only do sends if I have particles
if (nlocal_ > 0)
{
@ -185,29 +220,62 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
{
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
int tag = (me_ << 16) | dst_procs_[i];
#ifdef O2ODEBUG
std::cout<< "[" << me_ << "]"
<< " SEND to: " << dst_procs_[i]
<< " nlocal_ " << nlocal_
<< " data_length " << data_length
<< " tag " << tag
<< std::endl;
#endif
MPI_Isend
(
src,
nlocal_*data_length,
wrap.mpi_type,
dst_procs_[i],
0,
comm_
tag,
comm_,
&request_[requesti]
);
requesti++;
}
#ifdef O2ODEBUG
else
{
std::cout<< "[" << me_ << "]"
<< " SEND to: " << dst_procs_[i]
<< " skipped (self)" << std::endl;
}
#endif
}
}
#ifdef O2ODEBUG
else
{
std::cout<< "[" << me_ << "]" << " SEND skipped (no data)" << std::endl;
}
#endif
// only wait if requests were actually posted
if (requesti > 0)
{
#ifdef O2ODEBUG
std::cout<< "[" << me_ << "]"
<< " WAIT for: " << requesti << " requests" << std::endl;
#endif
MPI_Waitall(requesti, request_, status_);
}
#ifdef O2ODEBUG
else
{
std::cout<< "[" << me_ << "]" << " WAIT skipped (no requests)" << std::endl;
}
#endif
// copy on-proc data
if (offset_local > -1)
@ -223,6 +291,13 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
dst[locali+offset_local] = src[locali];
}
}
// wait for all procs to complete data exchange
#ifdef O2ODEBUG
std::cout<< "[" << me_ << "] BARRIER" << std::endl;
#endif
MPI_Barrier(comm_);
}
template void One2One::exchange<int>(int*&, int*&, int);

7
src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
View File

@ -183,6 +183,7 @@ void twoWayOne2One::createProcMap()
point(ligbb[0][0], ligbb[0][1], ligbb[0][2]),
point(ligbb[1][0], ligbb[1][1], ligbb[1][2])
);
delete [] ligbb;
if (boundBoxMargin_ > 0.)
{
thisLigBox.max() = thisLigBox.max() + boundBoxMargin_ * vector::one;
@ -230,6 +231,7 @@ twoWayOne2One::~twoWayOne2One()
{
delete foam2lig_;
delete lig2foam_;
delete [] lig2foam_mask_;
destroy(lig2foam_ids_);
destroy(foam2lig_ids_);
@ -777,10 +779,7 @@ void twoWayOne2One::locateParticles()
lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids_);
}
if (lig2foam_mask_)
{
delete [] lig2foam_mask_;
}
delete [] lig2foam_mask_;
lig2foam_mask_ = new bool[lig2foam_->ncollected_];
DynamicList<label> cellIds;

3
src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H
View File

@ -105,8 +105,6 @@ public:
virtual void postFlow() {}
virtual void solve() {}
const volScalarField& kf0Field() const;
const volScalarField& CpField() const;
@ -116,7 +114,6 @@ public:
return 0.0;
}
};

4
src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferInterGrain/heatTransferInterGrain.C
View File

@ -283,8 +283,6 @@ void heatTransferInterGrain::calcEnergyContribution()
if(cellI >= 0)
{
voidfraction = voidfraction_[cellI];
if (voidfraction < 0.01)
voidfraction = 0.01;
partVolume = particleCloud_.particleVolume(index);
QPartPart = QPartPart_[cellI];
@ -426,8 +424,6 @@ void heatTransferInterGrain::calcPartThermRad()
scalar heatTransferInterGrain::FE(scalar voidfraction, scalar emissivity, scalar L)
{
if (voidfraction < 1 - SMALL) return 1.0; // limit of voidfraction -> 1.0
scalar B = 1.25 * pow((1-voidfraction)/voidfraction,1.111);
scalar sqrt_aP = sqrt(1-voidfraction);
scalar x = 2.0/emissivity - 1.0;

2
src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferInterGrain/heatTransferInterGrain.H
View File

@ -104,7 +104,7 @@ protected:
const scalar TMin = 10.0;
const scalar voidfracMax = 0.95;
const scalar voidfracMax = 0.8; // maximum voidfraction for which B-B correlation is valid; TODO: check value!
const scalar voidfracMin = 0.05;

1
src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C
View File

@ -77,7 +77,6 @@ BeetstraDrag::BeetstraDrag
particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force
particleCloud_.probeM().vectorFields_.append("Urel");
particleCloud_.probeM().scalarFields_.append("Rep");
particleCloud_.probeM().scalarFields_.append("betaP");
particleCloud_.probeM().scalarFields_.append("voidfraction");
particleCloud_.probeM().writeHeader();

313
src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.C
View File

@ -34,7 +34,6 @@ Description
#include "ParmarBassetForce.H"
#include "addToRunTimeSelectionTable.H"
#include "smoothingModel.H"
#include "constDiffSmoothing.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -71,17 +70,18 @@ ParmarBassetForce::ParmarBassetForce
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
UsFieldName_(propsDict_.lookup("granVelFieldName")),
Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
nInt_(readLabel(propsDict_.lookup("nIntegral"))),
discOrder_(readLabel(propsDict_.lookup("discretisationOrder"))),
nHist_(nInt_+discOrder_+1),
FHistSize_(2*discOrder_+1),
ddtUrelHist_(nHist_,NULL), // UrelHist_[ndt in past][particle ID][dim]
rHist_(nHist_,NULL), // rHist_[ndt in past][particle ID][0]
FHist_(2,List<double**>(FHistSize_,NULL)), // FHist_[k={1,2}-1][ndt in past][particle ID][dim]
gH0_(NULL),
tRef_(NULL),
mRef_(NULL),
lRef_(NULL),
nInt_(propsDict_.lookupOrDefault<int>("nIntegral", -1)),
discOrder_(propsDict_.lookupOrDefault<int>("discretisationOrder", 1)),
ddtUrelHistSize_(nInt_+discOrder_),
rHistSize_(nInt_),
FHistSize_(2*discOrder_),
ddtUrelHistRegName_(typeName + "ddtUrelHist"), // indexed as: ddtUrelHist_[particle ID][iDim*ddtUrelHistSize_+iHist]
rHistRegName_(typeName + "rHist"), // indexed as: rHist_[particle ID][iHist]
FHistRegName_(typeName + "FHist"), // indexed as: ddtUrelHist_[particle ID][iDim*FHistSize_*2+iHist*2+k]
gH0RegName_(typeName + "gH0"),
tRefRegName_(typeName + "tRef"),
mRefRegName_(typeName + "mRef"),
lRefRegName_(typeName + "lRef"),
Urel_
( IOobject
(
@ -115,12 +115,12 @@ ParmarBassetForce::ParmarBassetForce
)
)
{
// init force sub model
setForceSubModels(propsDict_);
// define switches which can be read from dict
forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
forceSubM(0).readSwitches();
//Extra switches/settings
@ -129,15 +129,23 @@ ParmarBassetForce::ParmarBassetForce
if (discOrder_ < 1 || discOrder_ > 2)
FatalError << "Parmar Basset Force: Discretisation order > 2 not implemented!" << abort(FatalError);
if (nInt_ < 1)
FatalError << "Parmar Basset Force: nIntegral missing or invalid, must be > 1" << abort(FatalError);
// allocate particle properties
particleCloud_.registerParticleProperty<double**>(ddtUrelHistRegName_,3*ddtUrelHistSize_,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( rHistRegName_, rHistSize_,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( FHistRegName_, 6*FHistSize_,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( gH0RegName_, 1 ,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( tRefRegName_, 1 ,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( mRefRegName_, 1 ,NOTONCPU,false);
particleCloud_.registerParticleProperty<double**>( lRefRegName_, 1 ,NOTONCPU,false);
//Append the field names to be probed
particleCloud_.probeM().initialize(typeName, typeName+".logDat");
particleCloud_.probeM().vectorFields_.append("ParmarBassetForce"); //first entry must the be the force
particleCloud_.probeM().vectorFields_.append("Urel");
particleCloud_.probeM().vectorFields_.append("ddtUrel");
//
particleCloud_.probeM().vectorFields_.append("UrelNoSmooth");
particleCloud_.probeM().vectorFields_.append("ddtUrelNoSmooth");
//
particleCloud_.probeM().vectorFields_.append("Fshort");
particleCloud_.probeM().vectorFields_.append("Flong1");
particleCloud_.probeM().vectorFields_.append("Flong2");
@ -158,20 +166,6 @@ ParmarBassetForce::ParmarBassetForce
ParmarBassetForce::~ParmarBassetForce()
{
particleCloud_.dataExchangeM().destroy(gH0_, 1);
particleCloud_.dataExchangeM().destroy(tRef_, 1);
particleCloud_.dataExchangeM().destroy(mRef_, 1);
particleCloud_.dataExchangeM().destroy(lRef_, 1);
for (int i=0; i<nHist_; i++)
{
particleCloud_.dataExchangeM().destroy(ddtUrelHist_[i],3);
particleCloud_.dataExchangeM().destroy(rHist_ [i],1);
}
for (int i=0; i<FHistSize_; i++)
for (int k=0; k<2; k++)
particleCloud_.dataExchangeM().destroy(FHist_[k][i],3);
}
@ -180,16 +174,18 @@ ParmarBassetForce::~ParmarBassetForce()
void ParmarBassetForce::setForce() const
{
// allocate arrays
if(particleCloud_.numberOfParticlesChanged())
reAllocArrays();
double**& ddtUrelHist_ = particleCloud_.getParticlePropertyRef<double**>(ddtUrelHistRegName_);
double**& rHist_ = particleCloud_.getParticlePropertyRef<double**>( rHistRegName_);
double**& FHist_ = particleCloud_.getParticlePropertyRef<double**>( FHistRegName_);
double**& gH0_ = particleCloud_.getParticlePropertyRef<double**>( gH0RegName_);
double**& tRef_ = particleCloud_.getParticlePropertyRef<double**>( tRefRegName_);
double**& mRef_ = particleCloud_.getParticlePropertyRef<double**>( mRefRegName_);
double**& lRef_ = particleCloud_.getParticlePropertyRef<double**>( lRefRegName_);
vector position(0,0,0);
vector Urel(0,0,0);
vector ddtUrel(0,0,0);
scalar t0min = 0.00;
const volScalarField& nufField = forceSubM(0).nuField();
const volScalarField& rhoField = forceSubM(0).rhoField();
@ -220,26 +216,17 @@ void ParmarBassetForce::setForce() const
Urel_ = Us_ - U_;
ddtUrel_ = fvc::ddt(Us_) - fvc::ddt(U_) - (Us_ & fvc::grad(U_));
//
volVectorField UrelNoSmooth_ = Urel_;
volVectorField ddtUrelNoSmooth_ = ddtUrel_;
//
smoothingM().smoothen(Urel_);
smoothingM().smoothen(ddtUrel_);
interpolationCellPoint<vector> UrelInterpolator_(Urel_);
interpolationCellPoint<vector> ddtUrelInterpolator_(ddtUrel_);
//
interpolationCellPoint<vector> UrelNoSmoothInterpolator_(UrelNoSmooth_);
interpolationCellPoint<vector> ddtUrelNoSmoothInterpolator_(ddtUrelNoSmooth_);
//
for(int index = 0;index < particleCloud_.numberOfParticles(); index++)
{
vector ParmarBassetForce(0,0,0);
vector Fshort(0,0,0);
vector Flong[2]={vector::zero, vector::zero};
label cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1) // particle Found
@ -268,12 +255,13 @@ void ParmarBassetForce::setForce() const
if (gH0_[index][0]!=NOTONCPU)
{
r = pow(gH0_[index][0]/gH,1.5); // Eq. 3.4
scalar gHratio = gH0_[index][0]/gH;
r = gHratio*sqrt(gHratio); // gHratio^1.5, Eq. 3.4
if (r<0.25 || r>2.0)
{
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
}
@ -285,7 +273,7 @@ void ParmarBassetForce::setForce() const
scalar Vs = rs*rs*rs*M_PI*4/3;
scalar mRef = Vs*rhoField[cellI] * gH * 5.2863; // 9/(2*sqrt(pi))*(256/pi)^(1/6) = 5.2863 (Eq. 3.2)
gH0_[index][0] = gH;
gH0_[index][0] = gH;
tRef_[index][0] = tRef;
mRef_[index][0] = mRef;
lRef_[index][0] = rs;
@ -298,29 +286,15 @@ void ParmarBassetForce::setForce() const
scalar dt = U_.mesh().time().deltaT().value() / tRef_[index][0]; // dim.less
scalar t0 = nInt_*dt; // dim.less
//********* update histories *********//
// non-dimensionlise
Urel /= mps;
ddtUrel /= mpss;
// update ddtUrel and r history
update_ddtUrelHist(ddtUrel,index); // add current dim.less ddtUrel to history
update_rHist(r,index); // add current r to history
// warning and reset for too small t0
if (t0<t0min)
{
Pout << "ParmarBassetForce WARNING: t0 = " << t0 << " at ID = " << index << endl;
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
// check length of known history
int nKnown = 0;
for (int j=0; j<nHist_; j++) // loop over past times
int nKnown = 1; // we always know the current step
for (int j=0; j<ddtUrelHistSize_; j++) // loop over past times
{
if (ddtUrelHist_[j][index][0] == NOTONCPU)
if (ddtUrelHist_[index][j] == NOTONCPU)
break;
else
nKnown++;
@ -332,117 +306,92 @@ void ParmarBassetForce::setForce() const
int nShort = min(nKnown,nInt_+1);
// int_0^dt K(r,xi) dxi * ddtU(t) dxi (singularity treated by assuming constant acceleration)
if (nShort>0)
{
for (int i=0; i<3; i++) // loop over dimensions
Fshort[i] = -calculateK0(r,dt) * ddtUrelHist_[0][index][i];
}
for (int i=0; i<3; i++) // loop over dimensions
Fshort[i] = -calculateK0(r,dt) * ddtUrel[i];
// int_dt^t0 K(r,xi) * ddtU(t-xi) dxi (trapezoid rule)
if (nShort>2)
{
for (int j=1; j<nShort; j++)
for (int j=0; j<(nShort-1); j++) // we don't use the current step here, hence nShort-1
{
scalar xi = j*dt;
scalar K = pow((pow(xi,.25) + rHist_[j][index][0]*xi),-2.); // Eq. 3.4
scalar xi = (j+1)*dt;
scalar invsqrtK = sqrt(sqrt(xi)) + rHist_[index][j]*xi; // K^-0.5
scalar K = 1./(invsqrtK*invsqrtK); // Eq. 3.4
for (int i=0; i<3; i++) // loop over dimensions
Fshort[i] -= trapWeight(j,nShort) * K * ddtUrelHist_[j][index][i] * dt;
Fshort[i] -= trapWeight(j,nShort-1) * K * ddtUrelHist_[index][i*ddtUrelHistSize_+j] * dt;
}
}
//********* long term force computing (differential form) *********//
// update F1, F2 history
update_FHist(vector::zero,vector::zero,index);
// initialise ddtUrel(t0) and Flong(:) as 0 and r(t0) as 1
if (nKnown == nInt_)
{
for (int j=nInt_; j<nHist_; j++) // loop over past times
{
rHist_[j][index][0] = 1.;
// initialise the histories beyond nInt
for (int j=nInt_-1; j<rHistSize_; j++) // loop over past times
rHist_[index][j] = 1.;
for (int j=nInt_-1; j<ddtUrelHistSize_; j++) // loop over past times
for (int i=0; i<3; i++) // loop over dimensions
ddtUrelHist_[j][index][i] = 0.0;
}
ddtUrelHist_[index][i*ddtUrelHistSize_+j] = 0.0;
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=0; j<FHistSize_; j++) // loop over past times
for (int i=0; i<3; i++) // loop over dimensions
FHist_[k][j][index][i] = 0.0;
nKnown = nHist_;
FHist_[index][i*FHistSize_*2+j*2+k] = 0.0;
nKnown = ddtUrelHistSize_+1;
}
// solve ODEs
if (nKnown == nHist_)
if (nKnown == ddtUrelHistSize_+1)
{
for (int k=0; k<2; k++) // loop over F1, F2
{
//calculate coefficients
double C[4];
calculateCoeffs(k,t0,rHist_[nInt_][index][0],c,chi,C);
calculateCoeffs(k,t0,rHist_[index][nInt_-1],c,chi,C);
// solve Eq. 3.20
solveFlongODE(k,C,dt,index);
Flong[k] = solveFlongODE(FHist_,ddtUrelHist_,k,C,dt,index);
}
}
//********* update histories *********//
update_ddtUrelHist(ddtUrelHist_,ddtUrel,index); // add current dim.less ddtUrel to history
update_rHist(rHist_,r,index); // add current r to history
update_FHist(FHist_,Flong[0],Flong[1],index);
//********* total force *********//
// sum and convert to N
for (int i=0; i<3; i++) // loop over dimensions
{
ParmarBassetForce[i] = Fshort[i];
for (int k=0; k<2; k++) // loop over F1, F2
ParmarBassetForce[i] += FHist_[k][0][index][i];
}
ParmarBassetForce = Fshort;
for (int k=0; k<2; k++) // loop over F1, F2
ParmarBassetForce += Flong[k];
ParmarBassetForce *= newton;
if (forceSubM(0).verbose() && index >= 0 && index < 2)
{
Pout << "cellI = " << cellI << endl;
Pout << "index = " << index << endl;
Pout << "Fshort = " << Fshort*newton << endl;
Pout << "Flong1 = " << Flong[0]*newton << endl;
Pout << "Flong2 = " << Flong[1]*newton << endl;
Pout << "Ftotal = " << ParmarBassetForce << endl;
}
// Set value fields and write the probe
if(probeIt_)
{
scalar ReRef = 0.75/(gH0_[index][0]-0.105);
vector Flong1(0,0,0);
vector Flong2(0,0,0);
for (int i=0; i<3; i++) // loop over dimensions
{
Flong1[i] = FHist_[0][0][index][i];
Flong2[i] = FHist_[1][0][index][i];
}
//
// relative velocity (m/s)
vector UrelNoSmooth;
vector ddtUrelNoSmooth;
if(forceSubM(0).interpolation())
UrelNoSmooth = UrelNoSmoothInterpolator_.interpolate(position,cellI);
else
UrelNoSmooth = UrelNoSmooth_[cellI];
// acceleration (m/s2)
if(forceSubM(0).interpolation())
ddtUrelNoSmooth = ddtUrelNoSmoothInterpolator_.interpolate(position,cellI);
else
ddtUrelNoSmooth = ddtUrelNoSmooth_[cellI];
UrelNoSmooth /= mps;
ddtUrelNoSmooth /= mpss;
//
#include "setupProbeModelfields.H"
vValues.append(ParmarBassetForce); //first entry must the be the force
vValues.append(Urel);
vValues.append(ddtUrel);
//
vValues.append(UrelNoSmooth);
vValues.append(ddtUrelNoSmooth);
//
vValues.append(Fshort);
vValues.append(Flong1);
vValues.append(Flong2);
vValues.append(Flong[0]);
vValues.append(Flong[1]);
sValues.append(ReRef);
sValues.append(tRef_[index][0]);
sValues.append(mRef_[index][0]);
@ -456,8 +405,8 @@ void ParmarBassetForce::setForce() const
else
{
// not on CPU
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
gH0_[index][0] = NOTONCPU; //reset reference
ddtUrelHist_[index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
}
// write particle based data to global array
@ -465,32 +414,11 @@ void ParmarBassetForce::setForce() const
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::reAllocArrays() const
{
Pout << "ParmarBassetForce::reAllocArrays..." << endl;
particleCloud_.dataExchangeM().allocateArray(gH0_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(tRef_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(mRef_, NOTONCPU,1);
particleCloud_.dataExchangeM().allocateArray(lRef_, NOTONCPU,1);
for (int i=0; i<nHist_; i++)
{
particleCloud_.dataExchangeM().allocateArray(ddtUrelHist_[i],NOTONCPU,3);
particleCloud_.dataExchangeM().allocateArray(rHist_ [i],NOTONCPU,1);
}
for (int i=0; i<2*discOrder_+1; i++)
for (int k=0; k<2; k++)
particleCloud_.dataExchangeM().allocateArray(FHist_[k][i],0.0,3);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
{
scalar cbrtr = cbrt(r); // cube root of r
scalar gamma = cbrtr*pow(dt,0.25);
scalar gamma = cbrtr*sqrt(sqrt(dt));
/*
scalar K0 = 2./(9.*pow(r,0.666)) *
@ -513,44 +441,44 @@ scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar Foam::ParmarBassetForce::trapWeight(int i, int n) const
{
if ( (i==1) || (i==(n-1)) )
if ( (i==0) || (i==(n-1)) )
return 0.5;
else
return 1.0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_ddtUrelHist(const vector& ddtUrel, int index) const
void Foam::ParmarBassetForce::update_ddtUrelHist(double**& ddtUrelHist_, const vector& ddtUrel, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
for (int j=nHist_-1; j>0; j--) // loop over past times
ddtUrelHist_[j][index][i] = ddtUrelHist_[j-1][index][i];
for (int j=ddtUrelHistSize_-1; j>0; j--) // loop over past times
ddtUrelHist_[index][i*ddtUrelHistSize_+j] = ddtUrelHist_[index][i*ddtUrelHistSize_+j-1];
ddtUrelHist_[0][index][i] = ddtUrel[i];
ddtUrelHist_[index][i*ddtUrelHistSize_] = ddtUrel[i];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_rHist(scalar r, int index) const
void Foam::ParmarBassetForce::update_rHist(double**& rHist_, scalar r, int index) const
{
for (int j=nHist_-1; j>0; j--) // loop over past times
rHist_[j][index][0] = rHist_[j-1][index][0];
for (int j=rHistSize_-1; j>0; j--) // loop over past times
rHist_[index][j] = rHist_[index][j-1];
rHist_[0][index][0] = r;
rHist_[index][0] = r;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::update_FHist(const vector& F1, const vector& F2, int index) const
void Foam::ParmarBassetForce::update_FHist(double**& FHist_, const vector& F1, const vector& F2, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=FHistSize_-1; j>0; j--) // loop over past times
FHist_[k][j][index][i] = FHist_[k][j-1][index][i];
FHist_[index][i*FHistSize_*2+j*2+k] = FHist_[index][i*FHistSize_*2+(j-1)*2+k];
FHist_[0][0][index][i] = F1[i];
FHist_[1][0][index][i] = F2[i];
FHist_[index][i*FHistSize_*2 ] = F1[i];
FHist_[index][i*FHistSize_*2+1] = F2[i];
}
}
@ -564,18 +492,20 @@ void Foam::ParmarBassetForce::calculateCoeffs(int k, scalar t0, scalar r, double
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int index) const
vector Foam::ParmarBassetForce::solveFlongODE(double**& FHist_, double**& ddtUrelHist_, int k, double C[4], scalar dt, int index) const
{
vector Flong = vector::zero;
if (discOrder_==1)
{
for (int i=0; i<3; i++) // loop over dimensions
{
FHist_[k][0][index][i] =
Flong[i] =
(
( C[1]/dt+2/(dt*dt)) * FHist_[k][1][index][i]
- ( 1/(dt*dt)) * FHist_[k][2][index][i]
- (C[2]+C[3]/dt ) * ddtUrelHist_[nInt_ ][index][i]
+ ( C[3]/dt ) * ddtUrelHist_[nInt_+1][index][i]
( C[1]/dt+2/(dt*dt)) * FHist_[index][i*FHistSize_*2 +k]
- ( 1/(dt*dt)) * FHist_[index][i*FHistSize_*2+2+k]
- (C[2]+C[3]/dt ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_-1]
+ ( C[3]/dt ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_ ]
) / (C[0]+C[1]/dt+1/(dt*dt)); // Eq. 3.20 using first order temporal discretisation
}
}
@ -583,39 +513,20 @@ void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int i
{
for (int i=0; i<3; i++) // loop over dimensions
{
FHist_[k][0][index][i] =
Flong[i] =
(
( 4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[k][1][index][i]
- ( C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[k][2][index][i]
+ ( 8/(4*dt*dt)) * FHist_[k][3][index][i]
- ( 1/(4*dt*dt)) * FHist_[k][4][index][i]
( 4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[index][i*FHistSize_*2 +k]
- ( C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+2+k]
+ ( 8/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+4+k]
- ( 1/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+6+k]
- (C[2] + 3*C[3]/(2*dt) ) * ddtUrelHist_[nInt_ ][index][i]
+ ( 4*C[3]/(2*dt) ) * ddtUrelHist_[nInt_+1][index][i]
- ( C[3]/(2*dt) ) * ddtUrelHist_[nInt_+2][index][i]
- (C[2] + 3*C[3]/(2*dt) ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_-1]
+ ( 4*C[3]/(2*dt) ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_ ]
- ( C[3]/(2*dt) ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_+1]
) / (C[0] + 3*C[1]/(2*dt) + 9/(4*dt*dt)); // Eq. 3.20 using second order temporal discretisation
}
}
}
void Foam::ParmarBassetForce::rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const
{
for (int i=0; i<3; i++) // loop over dimensions
{
// rescale ddtU history
for (int j=0; j<nHist_; j++) // loop over past times
if (ddtUrelHist_[j][index][i] != NOTONCPU)
ddtUrelHist_[j][index][i] *= lScale/(tScale*tScale);
// rescale F1, F2 history
for (int k=0; k<2; k++) // loop over F1, F2
for (int j=0; j<FHistSize_; j++) // loop over past times
FHist_[k][j][index][i] *= mScale*lScale/(tScale*tScale);
}
// rescale r history
for (int j=0; j<nHist_; j++) // loop over past times
if (rHist_[j][index][0] != NOTONCPU)
rHist_[j][index][0] /= pow(rScale,1.5);
return Flong;
}
} // End namespace Foam

41
src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.H
View File

@ -59,35 +59,37 @@ class ParmarBassetForce
private:
dictionary propsDict_;
word velFieldName_;
const word velFieldName_;
const volVectorField& U_;
word UsFieldName_;
const word UsFieldName_;
const volVectorField& Us_;
label nInt_; //no of timesteps to solve full integral
const int nInt_; //no of timesteps to solve full integral
label discOrder_; //ODE discretisation order
const int discOrder_; //ODE discretisation order
label nHist_; //no of timesteps to save ddtUrel history for
const int ddtUrelHistSize_; //no of timesteps to save ddtUrel history for
label FHistSize_;
const int rHistSize_; //no of timesteps to save r history for
mutable List<double**> ddtUrelHist_;
const int FHistSize_; //no of timesteps to save Flong history for
mutable List<double**> rHist_;
const word ddtUrelHistRegName_;
mutable List<List<double**>> FHist_;
const word rHistRegName_;
mutable double** gH0_;
const word FHistRegName_;
mutable double** tRef_;
const word gH0RegName_;
mutable double** mRef_;
const word tRefRegName_;
mutable double** lRef_;
const word mRefRegName_;
const word lRefRegName_;
mutable volVectorField Urel_;
@ -101,17 +103,15 @@ private:
scalar trapWeight(int i, int n) const;
void update_ddtUrelHist(const vector& ddtUrel, int index) const;
void update_ddtUrelHist(double**& ddtUrelHist_, const vector& ddtUrel, int index) const;
void update_rHist(scalar r, int index) const;
void update_rHist(double**& rHist_, scalar r, int index) const;
void update_FHist(const vector& F1, const vector& F2, int index) const;
void update_FHist(double**& FHist_, const vector& F1, const vector& F2, int index) const;
void calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const;
void solveFlongODE(int k, double C[4], scalar dt, int index) const;
void rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const;
vector solveFlongODE(double**& FHist_, double**& ddtUrelHist_, int k, double C[4], scalar dt, int index) const;
public:
@ -136,9 +136,6 @@ public:
// Member Functions
void setForce() const;
void reAllocArrays() const;
inline const smoothingModel& smoothingM() const
{
return smoothingModel_;

45
src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C
View File

@ -5,7 +5,8 @@
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright 2012-2015 DCS Computing GmbH, Linz
Copyright 2015- JKU Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling.
@ -66,10 +67,14 @@ ShirgaonkarIB::ShirgaonkarIB
verbose_(propsDict_.found("verbose")),
twoDimensional_(propsDict_.found("twoDimensional")),
depth_(1),
multisphere_(propsDict_.found("multisphere")), // drag for a multisphere particle
useTorque_(propsDict_.found("useTorque")),
velFieldName_(propsDict_.lookup("velFieldName")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
pressureFieldName_(propsDict_.lookup("pressureFieldName")),
p_(sm.mesh().lookupObject<volScalarField> (pressureFieldName_))
p_(sm.mesh().lookupObject<volScalarField> (pressureFieldName_)),
solidVolFractionName_(multisphere_?propsDict_.lookup("solidVolFractionName"):word::null),
lambda_(multisphere_?sm.mesh().lookupObject<volScalarField>(solidVolFractionName_):volScalarField::null())
{
//Append the field names to be probed
particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@ -109,6 +114,7 @@ void ShirgaonkarIB::setForce() const
label cellI;
vector drag;
vector torque;
volVectorField h=forceSubM(0).IBDragPerV(U_,p_);
@ -119,6 +125,7 @@ void ShirgaonkarIB::setForce() const
//if(mask[index][0])
//{
drag=vector::zero;
torque=vector::zero;
for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
{
@ -128,6 +135,12 @@ void ShirgaonkarIB::setForce() const
if (cellI > -1) // particle Found
{
drag += h[cellI]*h.mesh().V()[cellI];
if(useTorque_)
{
vector rc = particleCloud_.mesh().C()[cellI];
vector positionCenter = particleCloud_.position(index);
torque += (rc - positionCenter)^h[cellI]*h.mesh().V()[cellI];
}
}
}
@ -150,10 +163,38 @@ void ShirgaonkarIB::setForce() const
<< "," << particleCloud_.impForces()[index][1]
<< "," << particleCloud_.impForces()[index][2]
<< endl;
if(useTorque_)
{
for(int j=0;j<3;j++) particleCloud_.DEMTorques()[index][j] = torque[j]; // Adding to the particle torque;
if(verbose_) Info << "DEMTorques = " << particleCloud_.DEMTorques()[index][0] << "," << particleCloud_.DEMTorques()[index][1] << "," << particleCloud_.DEMTorques()[index][2] << endl;
}
//}
}
// Calculate the force if the particle is multisphere template
if(multisphere_) calcForce();
}
void ShirgaonkarIB::calcForce() const
{
vector dragMS;
volVectorField h=forceSubM(0).IBDragPerV(U_,p_);
dragMS = vector::zero;
forAll(lambda_,cellI)
{
if(lambda_[cellI] > 0) dragMS += h[cellI]*h.mesh().V()[cellI];
else dragMS = dragMS;
}
Pout << "Drag force on particle clump = " << dragMS[0] << ", " << dragMS[1] << ", " << dragMS[2] << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

13
src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H
View File

@ -5,7 +5,8 @@
www.cfdem.com
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright 2012-2015 DCS Computing GmbH, Linz
Copyright 2015- JKU Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling.
@ -66,6 +67,10 @@ private:
bool depth_;
const bool multisphere_;
const bool useTorque_;
word velFieldName_;
const volVectorField& U_;
@ -74,6 +79,12 @@ private:
const volScalarField& p_;
word solidVolFractionName_;
const volScalarField& lambda_;
void calcForce() const;
public:
//- Runtime type information

13
src/lagrangian/cfdemParticle/subModels/forceModel/gradPForceSmooth/gradPForceSmooth.C
View File

@ -61,6 +61,7 @@ gradPForceSmooth::gradPForceSmooth
U_(sm.mesh().lookupObject<volVectorField> (velocityFieldName_)),
useRho_(false),
useU_(false),
temporalSmoothing_(false),
addedMassCoeff_(0.0),
smoothingModel_
(
@ -96,7 +97,7 @@ gradPForceSmooth::gradPForceSmooth
if (modelType_ == "B")
{
FatalError <<"using model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
FatalError <<"using model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
}
else if (modelType_ == "Bfull")
{
@ -135,6 +136,11 @@ gradPForceSmooth::gradPForceSmooth
particleCloud_.probeM().scalarFields_.append("Vs");
particleCloud_.probeM().scalarFields_.append("rho");
particleCloud_.probeM().writeHeader();
if ((smoothingM().type() == "temporalSmoothing") || (smoothingM().type() == "constDiffAndTemporalSmoothing"))
{
temporalSmoothing_ = true;
}
}
@ -150,6 +156,11 @@ gradPForceSmooth::~gradPForceSmooth()
void gradPForceSmooth::setForce() const
{
volVectorField gradPField = fvc::grad(p_);
// if temporal smoothing is used, let pSmooth evolve on its own - smoothing model itself looks up p
// else, set pSmooth to current p and apply spatial smoothing
if (!temporalSmoothing_) pSmooth_ = p_;
if (pFieldName_ == "p_rgh")
{
const volScalarField& rho_ = particleCloud_.mesh().lookupObject<volScalarField>("rho");

2
src/lagrangian/cfdemParticle/subModels/forceModel/gradPForceSmooth/gradPForceSmooth.H
View File

@ -62,6 +62,8 @@ private:
bool useU_; // if false: substitution p=0.5*rho*U^2
bool temporalSmoothing_;
mutable double addedMassCoeff_; //added mass coefficient
autoPtr<smoothingModel> smoothingModel_;

2
src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C
View File

@ -103,7 +103,6 @@ interface::~interface()
void interface::setForce() const
{
Info << "interface::setForce" << endl;
for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
{
//if(mask[index][0])
@ -193,7 +192,6 @@ Info << "interface::setForce" << endl;
} // end if particle found on proc domain
//}// end if in mask
}// end loop particles
Info << "interface::setForce - done" << endl;
}

129
src/lagrangian/cfdemParticle/subModels/forceModel/transferFluidProperties/transferFluidProperties.C Normal file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
License
This is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This code is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this code. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
transfer fluid properties to LIGGGHTS
SourceFiles
transferFluidProperties.C
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "transferFluidProperties.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(transferFluidProperties, 0);
addToRunTimeSelectionTable
(
forceModel,
transferFluidProperties,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
transferFluidProperties::transferFluidProperties
(
const dictionary& dict,
cfdemCloud& sm
)
:
forceModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props"))
{
particleCloud_.registerParticleProperty<double**>("fluidDensity",1);
particleCloud_.registerParticleProperty<double**>("fluidViscosity",1);
// init force sub model
setForceSubModels(propsDict_);
// define switches which can be read from dict
forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
forceSubM(0).readSwitches();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
transferFluidProperties::~transferFluidProperties()
{
}
// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * //
void transferFluidProperties::setForce() const
{
double**& fluidDensity_ = particleCloud_.getParticlePropertyRef<double**>("fluidDensity");
double**& fluidViscosity_ = particleCloud_.getParticlePropertyRef<double**>("fluidViscosity");
const volScalarField& rhoField = forceSubM(0).rhoField();
const volScalarField& nufField = forceSubM(0).nuField();
interpolationCellPoint<scalar> rhoInterpolator_(rhoField);
interpolationCellPoint<scalar> nufInterpolator_(nufField);
label cellI = 0;
double rho = 0.;
double nuf = 0.;
vector position(0,0,0);
for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
{
cellI = particleCloud_.cellIDs()[index][0];
if (cellI >= 0)
{
if(forceSubM(0).interpolation())
{
position = particleCloud_.position(index);
rho = rhoInterpolator_.interpolate(position,cellI);
nuf = nufInterpolator_.interpolate(position,cellI);
}
else
{
rho = rhoField[cellI];
nuf = nufField[cellI];
}
fluidDensity_[index][0] = rho;
fluidViscosity_[index][0] = nuf*rho;
}
}
particleCloud_.dataExchangeM().giveData("fluidDensity","scalar-atom",fluidDensity_);
particleCloud_.dataExchangeM().giveData("fluidViscosity","scalar-atom",fluidViscosity_);
if (forceSubM(0).verbose()) Info << "give data done" << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

79
src/lagrangian/cfdemParticle/subModels/forceModel/transferFluidProperties/transferFluidProperties.H Normal file
View File

@ -0,0 +1,79 @@
/*---------------------------------------------------------------------------*\
License
This is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This code is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this code. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Description
transfer fluid properties to LIGGGHTS
SourceFiles
transferFluidProperties.C
\*---------------------------------------------------------------------------*/
#ifndef transferFluidProperties_H
#define transferFluidProperties_H
#include "forceModel.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class transferFluidProperties Declaration
\*---------------------------------------------------------------------------*/
class transferFluidProperties
:
public forceModel
{
private:
dictionary propsDict_;
public:
//- Runtime type information
TypeName("transferFluidProperties");
// Constructors
//- Construct from components
transferFluidProperties
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~transferFluidProperties();
// Member Functions
void setForce() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

7
src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C
View File

@ -99,7 +99,12 @@ virtualMassForce::virtualMassForce
)
)
{
particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false);
// allocate UrelOld only if needed
int UrelOldSize = 0;
if(!splitUrelCalculation_ || !useUs_ )
UrelOldSize = 3;
particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,UrelOldSize,NOTONCPU,false);
// init force sub model
setForceSubModels(propsDict_);

2
src/lagrangian/cfdemParticle/subModels/massTransferModel/massTransferModel/massTransferModel.H
View File

@ -105,8 +105,6 @@ public:
virtual void postFlow() {}
virtual void solve() {}
const volScalarField& D0Field() const;
const volScalarField& CsField() const;

2
src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C
View File

@ -219,7 +219,7 @@ void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> v
*sPtr << setprecision(writePrecision_);
forAll(vValues, iter)
{
// if(!probeDebug_ && iter>0) break;
if(!probeDebug_ && iter>0) break;
*sPtr << vValues[iter][0] << " ";
*sPtr << vValues[iter][1] << " ";
*sPtr << vValues[iter][2] << " ";

97
src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C
View File

@ -68,6 +68,7 @@ constDiffAndTemporalSmoothing::constDiffAndTemporalSmoothing
smoothingLength_(dimensionedScalar("smoothingLength",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
refFieldName_(propsDict_.lookupOrDefault<word>("refField","")),
gamma_(readScalar(propsDict_.lookup("smoothingStrength"))),
correctBoundary_(propsDict_.lookupOrDefault<bool>("correctBoundary",false)),
verbose_(false)
@ -111,12 +112,43 @@ void constDiffAndTemporalSmoothing::smoothen(volScalarField& fieldSrc) const
double deltaT = sSmoothField.mesh().time().deltaTValue();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
// do smoothing
solve
(
fvm::ddt(sSmoothField)
-fvm::laplacian(DT_, sSmoothField)
);
if(refFieldName_.empty())
{
// do spatial smoothing
solve
(
fvm::ddt(sSmoothField)
-fvm::laplacian(DT_, sSmoothField)
);
// create fields for temporal smoothing
volScalarField refField = sSmoothField;
sSmoothField.ref()=fieldSrc.oldTime().internalField();
sSmoothField.correctBoundaryConditions();
sSmoothField.oldTime()=fieldSrc.oldTime();
sSmoothField.oldTime().correctBoundaryConditions();
// do temporal smoothing
solve
(
fvm::ddt(sSmoothField)
-
gamma_/deltaT * (refField - fvm::Sp(1.0,sSmoothField))
);
}
else
{
const volScalarField& refField = particleCloud_.mesh().lookupObject<volScalarField>(refFieldName_);
// do spatial and temporal smoothing
solve
(
fvm::ddt(sSmoothField)
-fvm::laplacian(DT_, sSmoothField)
-gamma_/deltaT * (refField - fvm::Sp(1.0,sSmoothField))
);
}
// bound sSmoothField_
forAll(sSmoothField,cellI)
@ -152,28 +184,43 @@ void constDiffAndTemporalSmoothing::smoothen(volVectorField& fieldSrc) const
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT.value();
// do spacial smoothing
solve
(
fvm::ddt(vSmoothField)
-fvm::laplacian(DT_, vSmoothField)
);
if(refFieldName_.empty())
{
// do spatial smoothing
solve
(
fvm::ddt(vSmoothField)
-fvm::laplacian(DT_, vSmoothField)
);
// create fields for temporal smoothing
volVectorField refField = vSmoothField;
// create fields for temporal smoothing
volVectorField refField = vSmoothField;
vSmoothField.ref()=fieldSrc.oldTime().internalField();
vSmoothField.correctBoundaryConditions();
vSmoothField.oldTime()=fieldSrc.oldTime();
vSmoothField.oldTime().correctBoundaryConditions();
vSmoothField.ref()=fieldSrc.oldTime().internalField();
vSmoothField.correctBoundaryConditions();
vSmoothField.oldTime()=fieldSrc.oldTime();
vSmoothField.oldTime().correctBoundaryConditions();
// do temporal smoothing
solve
(
fvm::ddt(vSmoothField)
-
gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
);
// do temporal smoothing
solve
(
fvm::ddt(vSmoothField)
-
gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
);
}
else
{
const volVectorField& refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
// do spatial and temporal smoothing
solve
(
fvm::ddt(vSmoothField)
-fvm::laplacian(DT_, vSmoothField)
-gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
);
}
// create temporally smoothened boundary field
if (correctBoundary_)

5
src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.H
View File

@ -6,7 +6,7 @@
Christoph Goniva, christoph.goniva@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright (C) 2013- Graz University of
Copyright (C) 2013- Graz University of
Technology, IPPT
-------------------------------------------------------------------------------
License
@ -65,6 +65,7 @@ private:
dimensionedScalar smoothingLength_;
dimensionedScalar smoothingLengthReferenceField_;
mutable dimensionedScalar DT_;
word refFieldName_;
scalar gamma_;
bool correctBoundary_;
bool verbose_;
@ -93,7 +94,7 @@ public:
bool doSmoothing() const;
//void dSmoothing(volScalarField&) const;
void smoothen(volScalarField&) const;
void smoothen(volVectorField&) const;

2
src/lagrangian/cfdemParticle/subModels/smoothingModel/temporalSmoothing/temporalSmoothing.C
View File

@ -130,7 +130,7 @@ void Foam::temporalSmoothing::smoothen(volVectorField& fieldSrc) const
vSmoothField.oldTime()=fieldSrc;
vSmoothField.oldTime().correctBoundaryConditions();
volVectorField refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
const volVectorField& refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
solve

1
src/lagrangian/cfdemParticleComp/Make/files
View File

@ -90,6 +90,7 @@ $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C
$(forceModels)/BeetstraDrag/BeetstraDrag.C
$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
$(forceModels)/dSauter/dSauter.C
$(forceModels)/transferFluidProperties/transferFluidProperties.C
$(forceModels)/Fines/Fines.C
$(forceModels)/Fines/FinesFields.C
$(forceModels)/Fines/FanningDynFines.C

19
src/recurrence/recModel/standardRecModel/standardRecModel.C
View File

@ -159,7 +159,24 @@ standardRecModel::standardRecModel
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
standardRecModel::~standardRecModel()
{}
{
const objectRegistry& objReg = base_.mesh().thisDb();
for(int j=0; j<volScalarFieldNames_.size(); j++)
{
objReg.checkOut(volScalarFieldList_[j][virtualTimeIndex]);
}
for(int j=0; j<volVectorFieldNames_.size(); j++)
{
objReg.checkOut(volVectorFieldList_[j][virtualTimeIndex]);
}
for(int j=0; j<surfaceScalarFieldNames_.size(); j++)
{
objReg.checkOut(surfaceScalarFieldList_[j][virtualTimeIndex]);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

27
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allclean.sh Executable file
View File

@ -0,0 +1,27 @@
#!/bin/bash
#------------------------------------------------------------------------------
# Allclean script for falling sphere test case
# clean CFD and DEM part
# Achuth N. Balachandran Nair - October 2018
#------------------------------------------------------------------------------
#- source environment variables
. ~/.bashrc
#- source CFDEM functions
source $CFDEM_PROJECT_DIR/etc/functions.sh
#- define variables
casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
#- clean up case
echo "deleting data at: $casePath"
source $WM_PROJECT_DIR/bin/tools/CleanFunctions
cd $casePath/CFD
cleanCase
rm -r $casePath/CFD/clockData
rm $casePath/DEM/post/*.*
#- preserve post directory
touch $casePath/DEM/post/.gitignore
echo "done"

25
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allrun.sh Executable file
View File

@ -0,0 +1,25 @@
#!/bin/bash
#------------------------------------------------------------------------------
# allrun script for falling sphere testcase
# run falling_sphere_two_way_coupling
# Achuth N. Balachandran Nair - Oct. 2018
#------------------------------------------------------------------------------
source $CFDEM_PROJECT_DIR/etc/functions.sh
#- define variables
casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
echo $casePath
# check if mesh was built
if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then
echo "mesh was built before - using old mesh"
else
echo "mesh needs to be built"
cd CFD/
blockMesh
fi
#- run parallel CFD-DEM
bash $casePath/parCFDDEMrun.sh

41
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/U Executable file
View File

@ -0,0 +1,41 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value $internalField;
}
channelWall
{
type zeroGradient;
}
outlet
{
type fixedValue;
value $internalField;
}
}
// ************************************************************************* //

46
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/Us Executable file
View File

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

36
tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/dSmoothing → tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/p Normal file → Executable file
View File

@ -1,38 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object dSmoothing;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 0 0 0 0 0];
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
side-walls
inlet
{
type fixedValue;
value $internalField;
}
channelWall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type zeroGradient;
type fixedValue;
value $internalField;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

28
tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/nut → tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/phiIB Normal file → Executable file
View File

@ -1,17 +1,16 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
object phiIB;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -21,24 +20,27 @@ internalField uniform 0;
boundaryField
{
side-walls
inlet
{
type fixedValue;
value uniform 0;
}
inlet
channelWall
{
type calculated;
value uniform 0;
type zeroGradient;
}
outlet
{
type calculated;
type fixedValue;
value uniform 0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

17
tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/CFD/0/rho → tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho
View File

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -16,16 +16,16 @@ FoamFile
dimensions [1 -3 0 0 0 0 0];
internalField uniform 1.58;
internalField uniform 1;
boundaryField
{
side-walls
inlet
{
type zeroGradient;
}
inlet
channelWall
{
type zeroGradient;
}
@ -36,5 +36,4 @@ boundaryField
}
}
// ************************************************************************* //

44
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/voidfraction Executable file
View File

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object voidfraction;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 1;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
channelWall
{
type zeroGradient;
}
frontAndBack
{
type zeroGradient;
}
}
// ************************************************************************* //

101
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/couplingProperties Normal file
View File

@ -0,0 +1,101 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object couplingProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// sub-models & settings
modelType none;
checkPeriodicCells;
verbose;
couplingInterval 10;
depth 0;
voidFractionModel IB;
locateModel engineIB;
meshMotionModel noMeshMotion;
regionModel allRegion;
dataExchangeModel twoWayMPI;
IOModel basicIO;
probeModel off;
averagingModel dilute;
clockModel off;
smoothingModel off;
forceModels
(
ShirgaonkarIB
ArchimedesIB
);
momCoupleModels
(
);
turbulenceModelType "turbulenceProperties";
// sub-model properties
ShirgaonkarIBProps
{
velFieldName "U";
pressureFieldName "p";
densityFieldName "rho";
verbose;
useTorque;
}
ArchimedesIBProps
{
densityFieldName "rho";
gravityFieldName "g";
voidfractionFieldName "voidfractionNext";
}
twoWayMPIProps
{
liggghtsPath "../DEM/in.liggghts_run";
}
IBProps
{
maxCellsPerParticle 20000;
alphaMin 0.30;
scaleUpVol 1.0;
}
engineIBProps
{
engineProps
{
treeSearch false;
}
zSplit 8;
xySplit 16;
}
// ************************************************************************* //

20
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/dynamicMeshDict Normal file
View File

@ -0,0 +1,20 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dynamicFvMesh staticFvMesh;
// ************************************************************************* //

21
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/g Executable file
View File

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

22
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/liggghtsCommands Executable file
View File

@ -0,0 +1,22 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object liggghtsCommands;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
liggghtsCommandModels
(
runLiggghts
);
// ************************************************************************* //

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