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 branch 'develop' of github.com:ParticulateFlow/CFDEMcoupling into develop

Browse Source
This commit is contained in:
tlichtenegger
2023-01-03 09:31:37 +01:00
parent 68318ab863 c05732bb45
commit ff42a04be1
34 changed files with 613 additions and 340 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/cfdemSolverMultiphaseScalar/Allwclean
View File

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

2
applications/solvers/cfdemSolverMultiphaseScalar/Allwmake
View File

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

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

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

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

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

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

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

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

@ -1,5 +1,5 @@
phase/phase.C phase/phase.C
alphaContactAngle/alphaContactAngleFvPatchScalarField.C alphaContactAngle/alphaContactAngleFvPatchScalarField.C
multiphaseMixture.C multiphaseMixtureScalar.C
LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoamScalar 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 Description
Contact-angle boundary condition for multi-phase interface-capturing Contact-angle boundary condition for multi-phase interface-capturing
simulations. Used in conjuction with multiphaseMixture. simulations. Used in conjuction with multiphaseMixtureScalar.
SourceFiles SourceFiles
alphaContactAngleFvPatchScalarField.C alphaContactAngleFvPatchScalarField.C
@ -37,7 +37,7 @@ SourceFiles
#define alphaContactAngleFvPatchScalarField_H #define alphaContactAngleFvPatchScalarField_H
#include "zeroGradientFvPatchFields.H" #include "zeroGradientFvPatchFields.H"
#include "multiphaseMixture.H" #include "multiphaseMixtureScalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -117,8 +117,8 @@ public:
typedef HashTable typedef HashTable
< <
interfaceThetaProps, interfaceThetaProps,
multiphaseMixture::interfacePair, multiphaseMixtureScalar::interfacePair,
multiphaseMixture::interfacePair::hash multiphaseMixtureScalar::interfacePair::hash
> thetaPropsTable; > 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 "alphaContactAngleFvPatchScalarField.H"
#include "Time.H" #include "Time.H"
#include "subCycle.H" #include "subCycle.H"
@ -31,13 +31,13 @@ License
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
const Foam::scalar Foam::multiphaseMixture::convertToRad = const Foam::scalar Foam::multiphaseMixtureScalar::convertToRad =
Foam::constant::mathematical::pi/180.0; Foam::constant::mathematical::pi/180.0;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::multiphaseMixture::calcAlphas() void Foam::multiphaseMixtureScalar::calcAlphas()
{ {
scalar level = 0.0; scalar level = 0.0;
alphas_ == 0.0; alphas_ == 0.0;
@ -51,7 +51,7 @@ void Foam::multiphaseMixture::calcAlphas()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::calcNu() const Foam::multiphaseMixtureScalar::calcNu() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -74,7 +74,7 @@ Foam::multiphaseMixture::calcNu() const
} }
Foam::tmp<Foam::surfaceScalarField> Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::calcStf() const Foam::multiphaseMixtureScalar::calcStf() const
{ {
tmp<surfaceScalarField> tstf tmp<surfaceScalarField> tstf
( (
@ -134,7 +134,7 @@ Foam::multiphaseMixture::calcStf() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::multiphaseMixture::multiphaseMixture Foam::multiphaseMixtureScalar::multiphaseMixtureScalar
( (
const volVectorField& U, const volVectorField& U,
const surfaceScalarField& phi, const surfaceScalarField& phi,
@ -230,7 +230,7 @@ Foam::multiphaseMixture::multiphaseMixture
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::rho() const Foam::multiphaseMixtureScalar::rho() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -247,7 +247,7 @@ Foam::multiphaseMixture::rho() const
Foam::tmp<Foam::scalarField> 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(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -264,9 +264,9 @@ Foam::multiphaseMixture::rho(const label patchi) const
Foam::tmp<Foam::volScalarField> 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(); return rho()*nu();
// PtrDictionary<phase>::const_iterator iter = phases_.begin(); // PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -283,7 +283,7 @@ Foam::multiphaseMixture::mu() const
Foam::tmp<Foam::scalarField> 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(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -306,7 +306,7 @@ Foam::multiphaseMixture::mu(const label patchi) const
Foam::tmp<Foam::surfaceScalarField> Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::muf() const Foam::multiphaseMixtureScalar::muf() const
{ {
return nuf()*fvc::interpolate(rho()); return nuf()*fvc::interpolate(rho());
@ -327,13 +327,13 @@ Foam::multiphaseMixture::muf() const
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nu() const Foam::multiphaseMixtureScalar::nu() const
{ {
return nu_; return nu_;
} }
Foam::tmp<Foam::scalarField> Foam::tmp<Foam::scalarField>
Foam::multiphaseMixture::nu(const label patchi) const Foam::multiphaseMixtureScalar::nu(const label patchi) const
{ {
//return nu_.boundaryField()[patchi]; //return nu_.boundaryField()[patchi];
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -355,7 +355,7 @@ Foam::multiphaseMixture::nu(const label patchi) const
Foam::tmp<Foam::surfaceScalarField> Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::nuf() const Foam::multiphaseMixtureScalar::nuf() const
{ {
//return muf()/fvc::interpolate(rho()); //return muf()/fvc::interpolate(rho());
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -374,7 +374,7 @@ Foam::multiphaseMixture::nuf() const
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cp() const Foam::multiphaseMixtureScalar::Cp() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -395,7 +395,7 @@ Foam::multiphaseMixture::Cp() const
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::kf() const Foam::multiphaseMixtureScalar::kf() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -417,7 +417,7 @@ Foam::multiphaseMixture::kf() const
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::D() const Foam::multiphaseMixtureScalar::D() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -439,7 +439,7 @@ Foam::multiphaseMixture::D() const
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::Cs() const Foam::multiphaseMixtureScalar::Cs() const
{ {
PtrDictionary<phase>::const_iterator iter = phases_.begin(); PtrDictionary<phase>::const_iterator iter = phases_.begin();
@ -456,7 +456,7 @@ Foam::multiphaseMixture::Cs() const
} }
Foam::tmp<Foam::surfaceScalarField> Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::diffusionCorrection() const Foam::multiphaseMixtureScalar::diffusionCorrection() const
{ {
surfaceScalarField numerator surfaceScalarField numerator
@ -517,7 +517,7 @@ Foam::multiphaseMixture::diffusionCorrection() const
return correction; return correction;
} }
void Foam::multiphaseMixture::solve() void Foam::multiphaseMixtureScalar::solve()
{ {
correct(); correct();
@ -570,7 +570,7 @@ void Foam::multiphaseMixture::solve()
} }
void Foam::multiphaseMixture::correct() void Foam::multiphaseMixtureScalar::correct()
{ {
forAllIter(PtrDictionary<phase>, phases_, iter) 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& alpha1,
const volScalarField& alpha2 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& alpha1,
const volScalarField& alpha2 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 // The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall. // velocity on the direction of the interface, parallel to the wall.
void Foam::multiphaseMixture::correctContactAngle void Foam::multiphaseMixtureScalar::correctContactAngle
( (
const phase& alpha1, const phase& alpha1,
const phase& alpha2, 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& alpha1,
const phase& alpha2 const phase& alpha2
@ -742,7 +742,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixture::nearInterface() const Foam::multiphaseMixtureScalar::nearInterface() const
{ {
tmp<volScalarField> tnearInt tmp<volScalarField> tnearInt
( (
@ -768,7 +768,7 @@ Foam::multiphaseMixture::nearInterface() const
} }
void Foam::multiphaseMixture::solveAlphas void Foam::multiphaseMixtureScalar::solveAlphas
( (
const scalar cAlpha const scalar cAlpha
) )
@ -901,7 +901,7 @@ void Foam::multiphaseMixture::solveAlphas
} }
bool Foam::multiphaseMixture::read() bool Foam::multiphaseMixtureScalar::read()
{ {
if (transportModel::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 Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
Class Class
multiphaseMixture multiphaseMixtureScalar
Description 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 which is an incompressible multi-phase mixture with built in solution
for the phase fractions with interface compression for interface-capturing. for the phase fractions with interface compression for interface-capturing.
It has been extended to include the void fraction in the volume fraction It has been extended to include the void fraction in the volume fraction
@ -33,11 +33,11 @@ Description
between each phase-pair. between each phase-pair.
SourceFiles SourceFiles
multiphaseMixture.C multiphaseMixtureScalar.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef multiphaseMixture_H #ifndef multiphaseMixtureScalar_H
#define multiphaseMixture_H #define multiphaseMixtureScalar_H
#include "incompressible/transportModel/transportModel.H" #include "incompressible/transportModel/transportModel.H"
#include "IOdictionary.H" #include "IOdictionary.H"
@ -52,10 +52,10 @@ namespace Foam
{ {
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class multiphaseMixture Declaration Class multiphaseMixtureScalar Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
class multiphaseMixture class multiphaseMixtureScalar
: :
public IOdictionary, public IOdictionary,
public transportModel public transportModel
@ -191,7 +191,7 @@ public:
// Constructors // Constructors
//- Construct from components //- Construct from components
multiphaseMixture multiphaseMixtureScalar
( (
const volVectorField& U, const volVectorField& U,
const surfaceScalarField& phi, const surfaceScalarField& phi,
@ -200,7 +200,7 @@ public:
//- Destructor //- 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 Description
Single incompressible phase derived from the phase-fraction. 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. simulations.
SourceFiles 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/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-I$(LIB_SRC)/sampling/lnInclude -I$(LIB_SRC)/sampling/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \ EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\ -L$(CFDEM_LIB_DIR)\

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)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-Wno-deprecated-copy
EXE_LIBS = \ EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\ -L$(CFDEM_LIB_DIR)\
@ -24,4 +25,4 @@ EXE_LIBS = \
-lfvOptions \ -lfvOptions \
-l$(CFDEM_LIB_NAME) \ -l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \ $(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS) $(CFDEM_ADD_LIBS)

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

@ -27,6 +27,7 @@ EXE_INC = \
-I$(LIB_SRC)/ODE/lnInclude \ -I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \ -I$(LIB_SRC)/combustionModels/lnInclude \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-Wno-deprecated-copy
EXE_LIBS = \ EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\ -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)/lagrangian/cfdemParticle/cfdTools \
-I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-Wno-deprecated-copy
EXE_LIBS = \ EXE_LIBS = \
-L$(CFDEM_LIB_DIR)\ -L$(CFDEM_LIB_DIR)\

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

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

10
doc/forceModel_ParmarBassetForce.txt
View File

@ -23,17 +23,19 @@ ParmarBassetForceProps
nIntegral scalar1; nIntegral scalar1;
discretisationOrder scalar2; discretisationOrder scalar2;
treatForceExplicit switch1; treatForceExplicit switch1;
interpolation switch2; verbose switch2;
interpolation switch3;
smoothingModel "smoothingModel"; smoothingModel "smoothingModel";
\} :pre \} :pre
{U} = name of the finite volume fluid velocity field :ulb,l {U} = name of the finite volume fluid velocity field :ulb,l
{Us} = name of the finite volume cell averaged particle velocity field :l {Us} = name of the finite volume cell averaged particle velocity field :l
{scalar1} = number of timesteps considered in the near history :l {scalar1} = number of timesteps considered in the near history (integer > 1):l
{scalar2} = (1 or 2) discretisation order of the far history differential equations :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 {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 {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 :ule
[Examples:] [Examples:]

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; smoothingModel constDiffSmoothing;
constDiffSmoothingProps constDiffSmoothingProps
\{ \{
lowerLimit number1; lowerLimit number1;
upperLimit number2; upperLimit number2;
smoothingLength lengthScale; smoothingLength lengthScale;
smoothingLengthReferenceField lengthScaleRefField; smoothingLengthReference lengthScaleRef;
smoothingLengthFieldName fieldName1;
smoothingLengthReferenceFieldName fieldName2;
verbose; verbose;
\} :pre \} :pre
{number1} = scalar fields will be bound to this lower value :ulb,l {number1} = scalar fields will be bound to this lower value :ulb,l
{number2} = scalar fields will be bound to this upper value :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 {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 {verbose} = (optional, default false) flag for debugging output :l
:ule :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 e.g. the average particle velocity, which are not specified in all cells in case
the flow is rather dilute. 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:] [Restrictions:]
This model is tested in a limited number of flow situations. This model is tested in a limited number of flow situations.

2
etc/library-list.txt
View File

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

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

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

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("dragForce"); //first entry must be the force
particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().vectorFields_.append("Urel");
particleCloud_.probeM().scalarFields_.append("Rep"); particleCloud_.probeM().scalarFields_.append("Rep");
particleCloud_.probeM().scalarFields_.append("betaP");
particleCloud_.probeM().scalarFields_.append("voidfraction"); particleCloud_.probeM().scalarFields_.append("voidfraction");
particleCloud_.probeM().writeHeader(); particleCloud_.probeM().writeHeader();

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

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

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

@ -59,35 +59,37 @@ class ParmarBassetForce
private: private:
dictionary propsDict_; dictionary propsDict_;
word velFieldName_; const word velFieldName_;
const volVectorField& U_; const volVectorField& U_;
word UsFieldName_; const word UsFieldName_;
const volVectorField& Us_; 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_; mutable volVectorField Urel_;
@ -101,17 +103,15 @@ private:
scalar trapWeight(int i, int n) const; 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 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; vector solveFlongODE(double**& FHist_, double**& ddtUrelHist_, int k, double C[4], scalar dt, int index) const;
void rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const;
public: public:
@ -136,9 +136,6 @@ public:
// Member Functions // Member Functions
void setForce() const; void setForce() const;
void reAllocArrays() const;
inline const smoothingModel& smoothingM() const inline const smoothingModel& smoothingM() const
{ {
return smoothingModel_; return smoothingModel_;

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 // init force sub model
setForceSubModels(propsDict_); setForceSubModels(propsDict_);

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_); *sPtr << setprecision(writePrecision_);
forAll(vValues, iter) forAll(vValues, iter)
{ {
// if(!probeDebug_ && iter>0) break; if(!probeDebug_ && iter>0) break;
*sPtr << vValues[iter][0] << " "; *sPtr << vValues[iter][0] << " ";
*sPtr << vValues[iter][1] << " "; *sPtr << vValues[iter][1] << " ";
*sPtr << vValues[iter][2] << " "; *sPtr << vValues[iter][2] << " ";

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

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

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

@ -159,7 +159,24 @@ standardRecModel::standardRecModel
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
standardRecModel::~standardRecModel() 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

8
tutorials/displacementField/CFD/constant/displacementProperties
View File

@ -31,6 +31,14 @@ dumpIndexDisplacementIncrement 8000;
dumpIndexInputIncrement 4000; dumpIndexInputIncrement 4000;
posIndex 0;
posX 1;
posY 2;
posZ 3;
nNeighMin 3; nNeighMin 3;
timePerDisplacementStep 0.1; timePerDisplacementStep 0.1;