diff --git a/README b/README new file mode 100755 index 00000000..fd2ca756 --- /dev/null +++ b/README @@ -0,0 +1,81 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright 2009-2012 JKU Linz + Copyright 2012- DCS Computing GmbH, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 3 of the License, or (at your + option) any later version. + + CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with CFDEMcoupling; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS + and OpenFOAM. Note: this code is not part of OpenFOAM (see DISCLAIMER). +\*---------------------------------------------------------------------------*/ + + +CFDEM coupling provides an open source parallel coupled CFD-DEM framework +combining the strengths of LIGGGHTS DEM code and the Open Source +CFD package OpenFOAM(R)(*). The CFDEMcoupling toolbox allows to expand +standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM +code LIGGGHTS. In this toolbox the particle representation within the +CFD solver is organized by "cloud" classes. Key functionalities are organised +in sub-models (e.g. force models, data exchange models, etc.) which can easily +be selected and combined by dictionary settings. + +The coupled solvers run fully parallel on distributed-memory clusters. + +Features are: + +- its modular approach allows users to easily implement new models +- its MPI parallelization enables to use it for large scale problems +- the "forum"_lws on CFD-DEM gives the possibility to exchange with other + users / developers +- the use of GIT allows to easily update to the latest version +- basic documentation is provided + +The file structure: + +- "src" directory including the source files of the coupling toolbox and models +- "applications" directory including the solver files for coupled CFD-DEM simulations +- "doc" directory including the documentation of CFDEMcoupling +- "tutorials" directory including basic tutorial cases showing the functionality + + + +Details on installation are given on the "www.cfdem.com" + +The functionality of this CFD-DEM framwork is described via "tutorial cases" showing +how to use different solvers and models. + +CFDEMcoupling stands for Computational Fluid Dynamics (CFD) - +Discrete Element Method (DEM) coupling. + +CFDEMcoupling is an open-source code, distributed freely under the terms of the +GNU Public License (GPL). + +Core development of CFDEMcoupling is done by +Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012 + + +\*---------------------------------------------------------------------------*/ +(*) "OpenFOAM(R)"_of is a registered trade mark of the ESI Group. +This offering is not affiliated, approved or endorsed by ESI Group, +the producer of the OpenFOAM® software and owner of the OpenFOAM® trade mark. +\*---------------------------------------------------------------------------*/ diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C index 76a6b58d..c07ca22d 100644 --- a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C +++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C @@ -83,95 +83,103 @@ int main(int argc, char *argv[]) particleCloud.clockM().stop("Coupling"); particleCloud.clockM().start(26,"Flow"); - // Pressure-velocity PISO corrector + + if(particleCloud.solveFlow()) { - // Momentum predictor - fvVectorMatrix UEqn - ( - fvm::ddt(voidfraction,U) + fvm::Sp(fvc::ddt(voidfraction),U) - + fvm::div(phi,U) + fvm::Sp(fvc::div(phi),U) -// + turbulence->divDevReff(U) - + particleCloud.divVoidfractionTau(U, voidfraction) - == - - fvm::Sp(Ksl/rho,U) - ); - - if (modelType=="B") - UEqn == - fvc::grad(p) + Ksl/rho*Us; - else - UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us; - - UEqn.relax(); - - if (momentumPredictor) - solve(UEqn); - - // --- PISO loop - - //for (int corr=0; corrdivDevReff(U) + + particleCloud.divVoidfractionTau(U, voidfraction) + == + - fvm::Sp(Ksl/rho,U) + ); if (modelType=="B") - U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA; + UEqn == - fvc::grad(p) + Ksl/rho*Us; else - U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA; + UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us; - U.correctBoundaryConditions(); + UEqn.relax(); - } // end piso loop + if (momentumPredictor) + solve(UEqn); + + // --- PISO loop + + //for (int corr=0; corrcorrect(); + }// end solveFlow + else + { + Info << "skipping flow solution." << endl; } - turbulence->correct(); - runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/doc/CFDEMcoupling_Manual.pdf b/doc/CFDEMcoupling_Manual.pdf index 2c69e8e2..56b4fead 100644 Binary files a/doc/CFDEMcoupling_Manual.pdf and b/doc/CFDEMcoupling_Manual.pdf differ diff --git a/doc/forceModel_GidaspowDrag.html b/doc/forceModel_GidaspowDrag.html index c31a57ce..0523427e 100644 --- a/doc/forceModel_GidaspowDrag.html +++ b/doc/forceModel_GidaspowDrag.html @@ -24,6 +24,7 @@ GidaspowDragProps voidfractionFieldName "voidfraction"; phi "scalar"; interpolation; + implDEM; };

Examples: diff --git a/doc/forceModel_GidaspowDrag.txt b/doc/forceModel_GidaspowDrag.txt index 5de896bf..83571706 100644 --- a/doc/forceModel_GidaspowDrag.txt +++ b/doc/forceModel_GidaspowDrag.txt @@ -22,6 +22,7 @@ GidaspowDragProps voidfractionFieldName "voidfraction"; phi "scalar"; interpolation; + implDEM; \}; :pre {U} = name of the finite volume fluid velocity field :ulb,l @@ -29,6 +30,7 @@ GidaspowDragProps {voidfraction} = name of the finite volume voidfraction field :l {phi} = drag correction factor (in doubt 1) :l {interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l +{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l :ule [Examples:] diff --git a/doc/forceModel_KochHillDrag.html b/doc/forceModel_KochHillDrag.html index 3938a881..db34b8aa 100644 --- a/doc/forceModel_KochHillDrag.html +++ b/doc/forceModel_KochHillDrag.html @@ -23,6 +23,7 @@ KochHillDragProps densityFieldName "density"; voidfractionFieldName "voidfraction"; interpolation; + implDEM; };

Examples: diff --git a/doc/forceModel_KochHillDrag.txt b/doc/forceModel_KochHillDrag.txt index d8ab7eb8..5a783fea 100644 --- a/doc/forceModel_KochHillDrag.txt +++ b/doc/forceModel_KochHillDrag.txt @@ -21,13 +21,14 @@ KochHillDragProps densityFieldName "density"; voidfractionFieldName "voidfraction"; interpolation; + implDEM; \}; :pre {U} = name of the finite volume fluid velocity field :ulb,l {density} = name of the finite volume gravity field :l {voidfraction} = name of the finite volume voidfraction field :l {interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l -{implDEM} = flag to use implicit formulation of drag on DEM side (normally off) :l +{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l :ule [Examples:] diff --git a/doc/githubAccess_public.pdf b/doc/githubAccess_public.pdf index 8de2281b..08d0fd52 100644 Binary files a/doc/githubAccess_public.pdf and b/doc/githubAccess_public.pdf differ diff --git a/doc/momCoupleModel.html b/doc/momCoupleModel.html index d40ba878..2acae6a2 100644 --- a/doc/momCoupleModel.html +++ b/doc/momCoupleModel.html @@ -29,6 +29,10 @@

Note: This examples list might not be complete - please look for other models (momCoupleModel_XY) in this documentation.

+

Forces can be coupled in an implicit way to the fluid solver (i.e., when solving the Navier-Stokes equations, the fluid velocity at the new time will be considered for the coupling force). This implicit coupling is typically done for the drag forces (look for "impForces()" in the implementation of the drag model). Implicit coupling is more stable (especially important for dense flows), but conflicts Newton's third law. Explicit forces are imposed on the flow solver in an explicit fashion (look for "expForces()" in the implementation of the drag model), which is less stable, but does not conflict Newton's third law. +

+

Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion. +

Description:

The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation. diff --git a/doc/momCoupleModel.txt b/doc/momCoupleModel.txt index cfbddc1b..05db72eb 100644 --- a/doc/momCoupleModel.txt +++ b/doc/momCoupleModel.txt @@ -27,6 +27,10 @@ momCoupleModels Note: This examples list might not be complete - please look for other models (momCoupleModel_XY) in this documentation. +Forces can be coupled in an implicit way to the fluid solver (i.e., when solving the Navier-Stokes equations, the fluid velocity at the new time will be considered for the coupling force). This implicit coupling is typically done for the drag forces (look for "impForces()" in the implementation of the drag model). Implicit coupling is more stable (especially important for dense flows), but conflicts Newton's third law. Explicit forces are imposed on the flow solver in an explicit fashion (look for "expForces()" in the implementation of the drag model), which is less stable, but does not conflict Newton's third law. + +Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion. + [Description:] The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation. diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 92a7d09f..83eff6da 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -6,7 +6,6 @@ voidFractionModels = subModels/voidFractionModel locateModels = subModels/locateModel meshMotionModels = subModels/meshMotionModel momCoupleModels = subModels/momCoupleModel -regionModels = subModels/regionModel dataExchangeModels = subModels/dataExchangeModel averagingModels = subModels/averagingModel clockModels = subModels/clockModel @@ -15,6 +14,7 @@ smoothingModels = subModels/smoothingModel probeModels = subModels/probeModel $(cfdemCloud)/cfdemCloud.C +derived/cfdemCloudBiDisperse/cfdemCloudBiDisperse.C derived/cfdemCloudIB/cfdemCloudIB.C derived/cfdemCloudMS/cfdemCloudMS.C @@ -23,36 +23,59 @@ $(forceModels)/forceModel/newForceModel.C $(forceModels)/noDrag/noDrag.C $(forceModels)/checkCouplingInterval/checkCouplingInterval.C $(forceModels)/DiFeliceDrag/DiFeliceDrag.C +$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C $(forceModels)/GidaspowDrag/GidaspowDrag.C $(forceModels)/SchillerNaumannDrag/SchillerNaumannDrag.C $(forceModels)/Archimedes/Archimedes.C $(forceModels)/ArchimedesIB/ArchimedesIB.C $(forceModels)/interface/interface.C $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C +$(forceModels)/interfaceParticleProbe/interfaceParticleProbe.C +$(forceModels)/fieldStore/fieldStore.C +$(forceModels)/fieldTimeAverage/fieldTimeAverage.C +$(forceModels)/fieldBound/fieldBound.C +$(forceModels)/volWeightedAverage/volWeightedAverage.C +$(forceModels)/totalMomentumExchange/totalMomentumExchange.C $(forceModels)/KochHillDrag/KochHillDrag.C $(forceModels)/KochHillRWDrag/KochHillRWDrag.C +$(forceModels)/BeetstraDrag/multiphaseFlowBasic/multiphaseFlowBasic.C +$(forceModels)/BeetstraDrag/BeetstraDrag.C +$(forceModels)/LaEuScalarLiquid/LaEuScalarLiquid.C $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C +$(forceModels)/LaEuScalarDust/LaEuScalarDust.C $(forceModels)/virtualMassForce/virtualMassForce.C $(forceModels)/gradPForce/gradPForce.C +$(forceModels)/gradULiftForce/gradULiftForce.C +$(forceModels)/HollowayDrag/HollowayDrag.C $(forceModels)/viscForce/viscForce.C $(forceModels)/MeiLift/MeiLift.C +$(forceModels)/melting/melting.C +$(forceModels)/KochHillDragNLift/KochHillDragNLift.C +$(forceModels)/stokesSpheroidDrag/stokesSpheroidDrag.C +$(forceModels)/solidsPressureForce/solidsPressureForce.C +$(forceModels)/periodicPressure/periodicPressure.C +$(forceModels)/periodicPressureControl/periodicPressureControl.C +$(forceModels)/averageSlipVel/averageSlipVel.C $(forceModels)/particleCellVolume/particleCellVolume.C -$(forceModels)/fieldTimeAverage/fieldTimeAverage.C -$(forceModels)/volWeightedAverage/volWeightedAverage.C $(forceModelsMS)/forceModelMS/forceModelMS.C $(forceModelsMS)/forceModelMS/newForceModelMS.C $(forceModelsMS)/DiFeliceDragMS/DiFeliceDragMS.C +$(forceModelsMS)/GidaspowDragMS/GidaspowDragMS.C +$(forceModelsMS)/noDragMS/noDragMS.C $(probeModels)/probeModel/probeModel.C $(probeModels)/probeModel/newProbeModel.C $(probeModels)/noProbe/noProbe.C $(probeModels)/particleProbe/particleProbe.C +/*$(probeModels)/interfaceParticleProbe/interfaceParticleProbe.C*/ $(IOModels)/IOModel/IOModel.C $(IOModels)/IOModel/newIOModel.C $(IOModels)/noIO/noIO.C $(IOModels)/basicIO/basicIO.C +$(IOModels)/tempIO/tempIO.C +$(IOModels)/colorIO/colorIO.C $(IOModels)/trackIO/trackIO.C $(IOModels)/sophIO/sophIO.C @@ -60,34 +83,36 @@ $(voidFractionModels)/voidFractionModel/voidFractionModel.C $(voidFractionModels)/voidFractionModel/newVoidFractionModel.C $(voidFractionModels)/centreVoidFraction/centreVoidFraction.C $(voidFractionModels)/dividedVoidFraction/dividedVoidFraction.C +$(voidFractionModels)/dividedVoidFractionBiDi/dividedVoidFractionBiDi.C $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C +$(voidFractionModels)/weightedNeigbhorsVoidFraction/weightedNeigbhorsVoidFraction.C $(locateModels)/locateModel/locateModel.C $(locateModels)/locateModel/newLocateModel.C $(locateModels)/standardSearch/standardSearch.C $(locateModels)/engineSearch/engineSearch.C +$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C $(locateModels)/turboEngineSearch/turboEngineSearch.C $(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C $(locateModels)/engineSearchIB/engineSearchIB.C - +$(locateModels)/hyperEngineSearch/hyperEngineSearch.C +$(locateModels)/ijkSearch/ijkSearch.C $(meshMotionModels)/meshMotionModel/meshMotionModel.C $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C $(meshMotionModels)/noMeshMotion/noMeshMotion.C +$(meshMotionModels)/DEMdrivenMeshMotion/DEMdrivenMeshMotion.C $(momCoupleModels)/momCoupleModel/momCoupleModel.C $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C $(momCoupleModels)/explicitCouple/explicitCouple.C +$(momCoupleModels)/explicitCoupleSource/explicitCoupleSource.C $(momCoupleModels)/implicitCouple/implicitCouple.C $(momCoupleModels)/noCouple/noCouple.C -$(regionModels)/regionModel/regionModel.C -$(regionModels)/regionModel/newRegionModel.C -$(regionModels)/allRegion/allRegion.C - $(dataExchangeModels)/dataExchangeModel/dataExchangeModel.C $(dataExchangeModels)/dataExchangeModel/newDataExchangeModel.C $(dataExchangeModels)/oneWayVTK/oneWayVTK.C @@ -95,11 +120,13 @@ $(dataExchangeModels)/twoWayFiles/twoWayFiles.C $(dataExchangeModels)/noDataExchange/noDataExchange.C $(dataExchangeModels)/twoWayMPI/twoWayMPI.C $(dataExchangeModels)/twoWayM2M/twoWayM2M.C +$(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C $(averagingModels)/averagingModel/averagingModel.C $(averagingModels)/averagingModel/newAveragingModel.C $(averagingModels)/dilute/dilute.C $(averagingModels)/dense/dense.C +$(averagingModels)/denseBiDi/denseBiDi.C $(clockModels)/clockModel/clockModel.C $(clockModels)/clockModel/newClockModel.C @@ -108,6 +135,7 @@ $(clockModels)/noClock/noClock.C $(liggghtsCommandModels)/liggghtsCommandModel/liggghtsCommandModel.C $(liggghtsCommandModels)/liggghtsCommandModel/newLiggghtsCommandModel.C +$(liggghtsCommandModels)/colorParticles/colorParticles.C $(liggghtsCommandModels)/execute/execute.C $(liggghtsCommandModels)/runLiggghts/runLiggghts.C $(liggghtsCommandModels)/writeLiggghts/writeLiggghts.C @@ -117,5 +145,6 @@ $(smoothingModels)/smoothingModel/smoothingModel.C $(smoothingModels)/smoothingModel/newSmoothingModel.C $(smoothingModels)/noSmoothing/noSmoothing.C $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C +$(smoothingModels)/localPSizeDiffSmoothing/localPSizeDiffSmoothing.C -LIB = $(FOAM_USER_LIBBIN)/lib$(CFDEM_LIB_NAME) +LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME) diff --git a/src/lagrangian/cfdemParticle/Make/options b/src/lagrangian/cfdemParticle/Make/options index 29b3e02c..3e648ff0 100644 --- a/src/lagrangian/cfdemParticle/Make/options +++ b/src/lagrangian/cfdemParticle/Make/options @@ -14,6 +14,7 @@ EXE_INC = \ -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \ -I$(CFDEM_LIGGGHTS_SRC_DIR) \ -I$(CFDEM_M2MLIB_PATH) \ + -I$(CFDEM_Many2ManyLIB_PATH) \ -I$(CFDEM_SRC_DIR)/cfdTools \ LIB_LIBS = \ @@ -28,4 +29,10 @@ LIB_LIBS = \ -L$(CFDEM_LIGGGHTS_SRC_DIR) \ -l$(CFDEM_LIGGGHTS_LIB_NAME) \ -L$(CFDEM_M2MLIB_PATH) \ - -lcouple + -lcouple \ + -L$(CFDEM_Many2ManyLIB_PATH) \ + -lcoupleMany2Many \ + +/* add -I$(CFDEM_POEMSLIB_PATH) \ to EXE_INC */ +/* -L$(CFDEM_POEMSLIB_PATH) \ */ +/* -lpoems */ diff --git a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H index b6507256..e721ad62 100755 --- a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H +++ b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H @@ -1,4 +1,4 @@ -word CFDEMversion="cfdem-2.6.3"; +word CFDEMversion="cfdem-2.6.4"; word compatibleLIGGGHTSversion="3.0.1"; word OFversion="2.2.x-commit-61b850bc107bdd60bbf1bf9a6417b9faf701d128"; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 8684785e..28358540 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -73,6 +73,7 @@ Foam::cfdemCloud::cfdemCloud IOobject::NO_WRITE ) ), + solveFlow_(true), verbose_(false), ignore_(false), modelType_(couplingProperties_.lookup("modelType")), @@ -95,8 +96,10 @@ Foam::cfdemCloud::cfdemCloud momCoupleModels_(couplingProperties_.lookup("momCoupleModels")), liggghtsCommandModelList_(liggghtsCommandDict_.lookup("liggghtsCommandModels")), turbulenceModelType_(couplingProperties_.lookup("turbulenceModelType")), + cg_(1.), cgOK_(true), impDEMdrag_(false), + imExSplitFactor_(1.0), useDDTvoidfraction_(false), ddtVoidfraction_ ( @@ -212,7 +215,10 @@ Foam::cfdemCloud::cfdemCloud #include "versionInfo.H" Info << "If BC are important, please provide volScalarFields -imp/expParticleForces-" << endl; - + if (couplingProperties_.found("solveFlow")) + solveFlow_=Switch(couplingProperties_.lookup("solveFlow")); + if (couplingProperties_.found("imExSplitFactor")) + imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor")); if (couplingProperties_.found("verbose")) verbose_=true; if (couplingProperties_.found("ignore")) ignore_=true; if (turbulenceModelType_=="LESProperties") @@ -259,7 +265,7 @@ Foam::cfdemCloud::cfdemCloud } dataExchangeM().setCG(); - if (!cgOK_ && forceM(0).cg() > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError); + if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError); } // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * // @@ -505,11 +511,47 @@ bool Foam::cfdemCloud::evolve if(verbose_) Info << "setvoidFraction done." << endl; clockM().stop("setvoidFraction"); - // set particles velocity field + // set average particles velocity field clockM().start(20,"setVectorAverage"); setVectorAverages(); - clockM().stop("setVectorAverage"); + + //Smoothen "next" fields + smoothingM().dSmoothing(); + smoothingM().smoothen(voidFractionM().voidFractionNext()); + smoothingM().smoothenReferenceField(averagingM().UsNext()); + + clockM().stop("setVectorAverage"); + } + + //============================================ + //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD + // IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED + // QUANTITIES AT THE GRID! + Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl; + clockM().start(24,"interpolateEulerFields"); + + // update voidFractionField + alpha = voidFractionM().voidFractionInterp(); + if(dataExchangeM().couplingStep() < 2) + { + alpha.oldTime() = alpha; // supress volume src + alpha.oldTime().correctBoundaryConditions(); + } + alpha.correctBoundaryConditions(); + + // calc ddt(voidfraction) + calcDdtVoidfraction(alpha); + + // update mean particle velocity Field + Us = averagingM().UsInterp(); + Us.correctBoundaryConditions(); + + clockM().stop("interpolateEulerFields"); + //============================================ + + if(doCouple) + { // set particles forces clockM().start(21,"setForce"); if(verbose_) Info << "- setForce(forces_)" << endl; @@ -530,32 +572,7 @@ bool Foam::cfdemCloud::evolve giveDEMdata(); clockM().stop("giveDEMdata"); }//end dataExchangeM().couple() - Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl; - clockM().start(24,"interpolateEulerFields"); - // update smoothing model - smoothingM().dSmoothing(); - - //============================================ - // update voidFractionField V1 - alpha = voidFractionM().voidFractionInterp(); - smoothingM().smoothen(alpha); - if(dataExchangeM().couplingStep() < 2) - { - alpha.oldTime() = alpha; // supress volume src - alpha.oldTime().correctBoundaryConditions(); - } - alpha.correctBoundaryConditions(); - - // calc ddt(voidfraction) - calcDdtVoidfraction(alpha); - - // update particle velocity Field - Us = averagingM().UsInterp(); - smoothingM().smoothenReferenceField(Us); - Us.correctBoundaryConditions(); - //============================================ - clockM().stop("interpolateEulerFields"); if(verbose_){ #include "debugInfo.H" diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H index 51e3f7f9..36c2bed5 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -88,6 +88,8 @@ protected: IOdictionary liggghtsCommandDict_; + Switch solveFlow_; + bool verbose_; bool ignore_; @@ -132,10 +134,14 @@ protected: const word turbulenceModelType_; + mutable scalar cg_; + bool cgOK_; bool impDEMdrag_; + mutable scalar imExSplitFactor_; + bool useDDTvoidfraction_; mutable volScalarField ddtVoidfraction_; @@ -235,10 +241,20 @@ public: label liggghtsCommandModelIndex(word); + inline void setCG(double) const; + + inline const scalar& cg() const; + + inline const bool& impDEMdrag() const; + + inline const scalar& imExSplitFactor() const; + inline const bool& ignore() const; inline const fvMesh& mesh() const; + inline bool solveFlow() const; + inline bool verbose() const; inline const IOdictionary& couplingProperties() const; @@ -281,6 +297,10 @@ public: virtual inline int ** particleTypes() const {return NULL;}; virtual label particleType(label index) const {return -1;}; + //access to the particles' orientation information + virtual inline double ** exArray() const {return NULL;}; + virtual vector ex(int) const {return Foam::vector(0,0,0);}; + inline int numberOfParticles() const; inline bool numberOfParticlesChanged() const; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 698ac718..170b5261 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -44,6 +44,27 @@ namespace Foam { // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline void cfdemCloud::setCG(double cg) const +{ + cg_ = cg; + Info << "cg is set to: " << cg_ << endl; +}; + +inline const bool& cfdemCloud::impDEMdrag() const +{ + return impDEMdrag_; +}; + +inline const scalar& cfdemCloud::imExSplitFactor() const +{ + return imExSplitFactor_; +}; + +inline const scalar& cfdemCloud::cg() const +{ + return cg_; +}; + inline const bool& cfdemCloud::ignore() const { return ignore_; @@ -54,6 +75,11 @@ inline const fvMesh& cfdemCloud::mesh() const return mesh_; } +inline bool cfdemCloud::solveFlow() const +{ + return bool(solveFlow_); +} + inline bool cfdemCloud::verbose() const { return verbose_; diff --git a/src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H b/src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H index aa7bbb96..4531873a 100644 --- a/src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H +++ b/src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H @@ -1,5 +1,5 @@ -#define version23 // currently being tested -//#define version22 // currently being used +//#define version23 // currently being tested +#define version22 // currently being used //#define version21 //#define version16ext //#define version15 @@ -8,7 +8,7 @@ //#define comp // if comp is on - you must use Make/options_comp! //define multi sphere -#define multisphere +//#define multisphere // features of 2.1 work also in 2.3 #if defined(version23) diff --git a/src/lagrangian/cfdemParticle/etc/bashrc b/src/lagrangian/cfdemParticle/etc/bashrc index 65703dd4..35aae2de 100755 --- a/src/lagrangian/cfdemParticle/etc/bashrc +++ b/src/lagrangian/cfdemParticle/etc/bashrc @@ -43,6 +43,10 @@ export CFDEM_LIGGGHTS_LIB_NAME=lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME #- CFDEM lib name export CFDEM_LIB_NAME=lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION +#- LMP Many2Many lib path and makefile +export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library +export CFDEM_Many2ManyLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME + #- LMP M2M lib path and makefile export CFDEM_M2MLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library export CFDEM_M2MLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME @@ -135,9 +139,6 @@ alias cfdemCompCFDEMsol='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil #- shortcut to compile CFDEMcoupling utilities alias cfdemCompCFDEMuti='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh' -#- shortcut to compile couple library -alias cfdemCompM2M='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileM2Mlib.sh' - #- shortcut to test basic tutorials alias cfdemTestTUT='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh' @@ -152,6 +153,10 @@ export -f cfdemLiggghts cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; } export -f cfdemLiggghtsPar +#- shortcut to open files including a pattern +cfdemGrep() { grep -rl "$1" ./* | xargs gedit; } +export -f cfdemGrep + # check if the run directory exists if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then : diff --git a/src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh b/src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh index d7ce0e0b..ac9f18d8 100755 --- a/src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh +++ b/src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh @@ -15,7 +15,6 @@ logDir="log" cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc mkdir -p $logDir - #================================================================================# # compile src #================================================================================# diff --git a/src/lagrangian/cfdemParticle/etc/cshrc b/src/lagrangian/cfdemParticle/etc/cshrc index 72683047..b5726625 100755 --- a/src/lagrangian/cfdemParticle/etc/cshrc +++ b/src/lagrangian/cfdemParticle/etc/cshrc @@ -43,6 +43,10 @@ setenv CFDEM_LIGGGHTS_LIB_NAME lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME #- CFDEM lib name setenv CFDEM_LIB_NAME lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION +#- LMP Many2Many lib path and makefile +setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library +setenv CFDEM_Many2ManyLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME + #- LMP M2M lib path setenv CFDEM_M2MLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library setenv CFDEM_M2MLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME @@ -133,9 +137,6 @@ alias cfdemCompCFDEMsol 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil #- shortcut to compile CFDEMcoupling utilities alias cfdemCompCFDEMuti 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh' -#- shortcut to compile couple library -alias cfdemCompM2M 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileM2Mlib.sh' - #- shortcut to test basic tutorials alias cfdemTestTUT 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh' diff --git a/src/lagrangian/cfdemParticle/etc/library-list.txt b/src/lagrangian/cfdemParticle/etc/library-list.txt index 27e798d4..e34c2541 100644 --- a/src/lagrangian/cfdemParticle/etc/library-list.txt +++ b/src/lagrangian/cfdemParticle/etc/library-list.txt @@ -1 +1,11 @@ lagrangian/cfdemParticle/dir + +#=====================================================' +#- RADL +fvOptions/dir +cylPorousMedia/dir + +#=====================================================' +#- other +finiteVolume/dir + diff --git a/src/lagrangian/cfdemParticle/etc/solver-list.txt b/src/lagrangian/cfdemParticle/etc/solver-list.txt index 925e9976..3d6d0341 100644 --- a/src/lagrangian/cfdemParticle/etc/solver-list.txt +++ b/src/lagrangian/cfdemParticle/etc/solver-list.txt @@ -1,4 +1,13 @@ -cfdemSolverPisoMS/dir cfdemSolverPiso/dir cfdemSolverIB/dir cfdemSolverPisoScalar/dir +cfdemSolverPimpleImEx/dir +cfdemSolverIBInterLubrication/dir +cfdemSolverIBScalar/dir +cfdemSolverInterDyM/dir +cfdemSolverInterDyMPC/dir +cfdemSolverBubble/dir +cfdemSolverPisoMS/dir +cfdemSolverPimpleDyM_22x/dir +cfdemSolverPimpleDyMMS_22x/dir +cfdemSolverPimpleDyMScalar_22x/dir diff --git a/src/lagrangian/cfdemParticle/etc/tutorial-list.txt b/src/lagrangian/cfdemParticle/etc/tutorial-list.txt index d893faf1..08ee103d 100644 --- a/src/lagrangian/cfdemParticle/etc/tutorial-list.txt +++ b/src/lagrangian/cfdemParticle/etc/tutorial-list.txt @@ -7,15 +7,51 @@ #===================================================================# cfdemSolverPiso/settlingTestMPI/dir - cfdemSolverPiso/ErgunTestMPI/dir - cfdemSolverPiso/ErgunTestMPI_cgs/dir - cfdemSolverPiso/ErgunTestMPI_restart/dir - cfdemSolverIB/twoSpheresGlowinskiMPI/dir - cfdemSolverPisoScalar/packedBedTemp/dir -cfdemSolverPiso/ErgunTestCG/dir +#===================================================================# +# RADL +cfdemSolverPimpleImEx/settlingTestMPI/dir +cfdemSolverPimpleImEx/ErgunTestMPI/dir +#cfdemSolverPimpleImEx/crossFlow/dir +#cfdemSolverIB/periodicCase/dir +#cfdemSolverIB/cfdemIBPeriodicCubicalBox_fullyPeriodic/dir +#cfdemSolverIBInterLubrication/twoCoatedParticlesRelMotion_smallTest/dir +#cfdemSolverIBScalar/cfdemIBPeriodicCubicalBoxScalar/dir + +#===================================================================# +# NesteJacobs +#Projects/Neste/cfdemSolverBubble/3pFBreactor/dir +#Projects/Neste/cfdemSolverInterDyM/3pFBreactor/dir + +#===================================================================# +# not in release: + +#cfdemSolverPiso/settlingTestBigParticleMPI/dir +#cfdemSolverPiso/ErgunTestCG/dir +#cfdemSolverPiso/ErgunTestM2M/dir +#cfdemSolverPiso/HopperEmptying/dir + +cfdemSolverPimpleDyM/ErgunTestMPI/dir + +#cfdemSolverPisoMS/settlingTestMPI/dir +#cfdemSolverPisoMS/ErgunTestMPI/dir + +#cfdemSolverInterDyM/twoPhaseSettlingTest/dir +#cfdemSolverInterDyM/ErgunTestMPI/dir +#cfdemSolverInterDyM/granularPiston/dir +#cfdemSolverInterDyM/sugarNcoffee/dir + +#cfdemSolverBubble/ErgunTestMPI_pureLiquid/dir + +#- these examples are already designed for 2.3.x +#cfdemSolverInterDyMPC/sugarNcoffee/dir +#cfdemSolverInterDyMPC/granularPiston/dir +#cfdemSolverInterDyMPC/meltingPot/dir + + + diff --git a/src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C b/src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C index 527a0c10..c2fd0e6d 100644 --- a/src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C +++ b/src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C @@ -45,8 +45,10 @@ defineRunTimeSelectionTable(IOModel, dictionary); // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void IOModel::dumpDEMdata() const -{} +int IOModel::dumpDEMdata() const +{ + return -1; +} fileName IOModel::createTimeDir(fileName path) const { @@ -64,6 +66,60 @@ fileName IOModel::createLagrangianDir(fileName path) const return cfdemCloudDirPath; } +fileName IOModel::buildFilePath(word dirName) const +{ + // create file structure + fileName path(""); + if(parOutput_) + { + path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud"); + mkDir(path,0777); + } else + { + path=fileName("."/dirName); + mkDir(path,0777); + mkDir(fileName(path/"constant"),0777); + OFstream* stubFile = new OFstream(fileName(path/"particles.foam")); + delete stubFile; + } + return path; +} + +void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className,word finaliser) const +{ + vector vec; + OFstream* fileStream = new OFstream(fileName(path/name)); + *fileStream << "FoamFile\n"; + *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object "<< name <<";}\n"; + *fileStream << nPProc <<"\n"; + + + if(type!="origProcId")*fileStream << "(\n"; + else if(type=="origProcId") + { + if(nPProc>0) *fileStream <<"{0}"<< "\n"; + else *fileStream <<"{}"<< "\n"; + } + + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + if (particleCloud_.cellIDs()[index][0] > -1) // particle Found + { + if (type=="scalar"){ + *fileStream << array[index][0] << " \n"; + }else if (type=="position" || type=="vector"){ + for(int i=0;i<3;i++) vec[i] = array[index][i]; + *fileStream <<"( "<< vec[0] <<" "<","0"); - streamDataToPath(lagPath_, particleCloud_.velocities(), "v","vector","vectorField",""); - streamDataToPath(lagPath_, particleCloud_.radii(), "r","scalar","scalarField",""); + streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","vector","Cloud","0"); + streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField",""); + streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField",""); } + return nPProc_; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Private Member Functions -fileName basicIO::buildFilePath(word dirName) const -{ - // create file structure - fileName path(""); - if(parOutput_) - { - path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud"); - mkDir(path,0777); - } else - { - path=fileName("."/dirName); - mkDir(path,0777); - mkDir(fileName(path/"constant"),0777); - OFstream* stubFile = new OFstream(fileName(path/"particles.foam")); - delete stubFile; - } - return path; -} -void basicIO::streamDataToPath(fileName path, double** array,word name,word type,word className,word finaliser) const -{ - vector vec; - OFstream* fileStream = new OFstream(fileName(path/name)); - *fileStream << "FoamFile\n"; - *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object "<< name <<";}\n"; - *fileStream << nPProc_ <<"\n"; - //*fileStream << "(\n"; - - if(type!="origProcId")*fileStream << "(\n"; - else if(type=="origProcId") - { - if(nPProc_>0) *fileStream <<"{0}"<< "\n"; - else *fileStream <<"{}"<< "\n"; - } - - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) - { - if (particleCloud_.cellIDs()[index][0] > -1) // particle Found - { - if (type=="scalar"){ - *fileStream << array[index][0] << " \n"; - }else if (type=="position" || type=="vector"){ - for(int i=0;i<3;i++) vec[i] = array[index][i]; - *fileStream <<"( "<< vec[0] <<" "<modify->find_fix_property(charName,"property/atom","scalar",0,0,"cfd coupling",false); @@ -320,7 +319,7 @@ void Foam::twoWayM2M::allocateArray for (int j = 0; j < width; j++) array[i][j] = initVal; } -void Foam::twoWayM2M::destroy(double** array,int len) const +void inline Foam::twoWayM2M::destroy(double** array,int len) const { // destroy array[i] first? lmp->memory->destroy(array); @@ -356,7 +355,7 @@ void Foam::twoWayM2M::allocateArray for (int j = 0; j < width; j++) array[i][j] = initVal; } -void Foam::twoWayM2M::destroy(int** array,int len) const +void inline Foam::twoWayM2M::destroy(int** array,int len) const { // destroy array[i] first? lmp->memory->destroy(array); @@ -370,7 +369,7 @@ void Foam::twoWayM2M::allocateArray(double*& array, double initVal, int length) for (int i = 0; i < len; i++) array[i] = initVal; } -void Foam::twoWayM2M::destroy(double* array) const +void inline Foam::twoWayM2M::destroy(double* array) const { lmp->memory->destroy(array); } @@ -383,7 +382,7 @@ void Foam::twoWayM2M::allocateArray(int*& array, int initVal, int length) const for (int i = 0; i < len; i++) array[i] = initVal; } -void Foam::twoWayM2M::destroy(int* array) const +void inline Foam::twoWayM2M::destroy(int* array) const { lmp->memory->destroy(array); } @@ -409,7 +408,7 @@ bool Foam::twoWayM2M::couple() const if(particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep())) { - const char* command = particleCloud_.liggghtsCommand()[i]().command(); + const char* command = particleCloud_.liggghtsCommand()[i]().command(0); Info << "Executing command: '"<< command <<"'"<< endl; lmp->input->one(command); } @@ -485,6 +484,8 @@ void Foam::twoWayM2M::syncIDs() const int* id_lammpsSync=NULL; double** pos_lammpsSync=NULL; + bool pos_lammps_alloc_flag=false; + bool id_lammps_alloc_flag=false; if(firstRun_ || particleLost_) { @@ -503,9 +504,11 @@ void Foam::twoWayM2M::syncIDs() const for (int j=0;j<3;j++) id_lammpsVec_[i*3+j] = id_lammps_[i]*3+j; } - - // get access to "x" + destroy(pos_lammps_,0); + pos_lammps_=NULL; pos_lammps_ = (double **) lammps_extract_atom(lmp,"x"); + pos_lammps_alloc_flag=false; + id_lammps_alloc_flag=false; } else { @@ -532,6 +535,8 @@ void Foam::twoWayM2M::syncIDs() const // map data according to last TS id_lammps_=NULL; allocateArray(id_lammps_,-1.,nlocal_foam_); + id_lammps_alloc_flag=true; + allocateArray(tmpI_,-1.,nlocal_foam_); lmp2foam_->exchange(id_lammpsSync, tmpI_); for(int i=0;iexchange(pos_lammpsSync ? pos_lammpsSync[0] : NULL, tmp_); allocateArray(pos_lammps_,0,3,nlocal_foam_); + pos_lammps_alloc_flag=true; for(int i=0;isetup(nlocal_foam_*3,id_foamVec_,nlocal_lammps_*3,id_lammpsVec_); foam2lmp_->setup(nlocal_foam_,id_foam_,nlocal_lammps_,id_lammpsSync); } + if (id_lammps_alloc_flag) destroy(id_lammps_); id_lammps_=NULL; // free pointer from LIG + destroy(id_lammpsSync); id_lammpsSync=NULL; // free pointer from LIG + if (pos_lammps_alloc_flag) destroy(pos_lammps_,0); pos_lammps_ = NULL; // free pointer from LIG particleCloud_.clockM().stop("setup_Comm"); @@ -629,7 +639,7 @@ void Foam::twoWayM2M::syncIDs() const //Info << "update communication - done." << endl; } -void Foam::twoWayM2M::locateParticle(int* id_lammpsSync) const +void Foam::twoWayM2M::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const { #if defined(version21) @@ -918,6 +928,9 @@ void Foam::twoWayM2M::locateParticle(int* id_lammpsSync) const } i++; } + Foam::dataExchangeModel::destroy(id_foam_nowhere_all); + id_foam_nowhere_all=NULL; + if (id_lammps_alloc_flag) destroy(id_lammps_); // make cpy id_lammps_=NULL; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/twoWayM2M.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/twoWayM2M.H index 0bb06191..4553a307 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/twoWayM2M.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/twoWayM2M.H @@ -188,22 +188,22 @@ public: // double ** void allocateArray(double**&, double, int, int) const; void allocateArray(double**&, double, int,const char* ="nparticles") const; - void destroy(double**,int) const; + void inline destroy(double**,int) const; //============ // int ** void allocateArray(int**&, int, int, int) const; void allocateArray(int**&, int, int,const char* ="nparticles") const; - void destroy(int**,int) const; + void inline destroy(int**,int) const; //============== //============== // double * void allocateArray(double*&, double, int) const; - void destroy(double*) const; + void inline destroy(double*) const; //============== // int * void allocateArray(int*&, int, int) const; - void destroy(int*) const; + void inline destroy(int*) const; //============== bool couple() const; @@ -212,9 +212,9 @@ public: int getNumberOfClumps() const; void syncIDs() const; - void locateParticle(int*) const; + void locateParticle(int*, bool) const; word myType() const{return typeName; }; - void setCG() const { particleCloud_.forceM(0).setCG(lmp->force->cg()); }; + void setCG() const { particleCloud_.setCG(lmp->force->cg()); }; }; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C index 8bc84fcf..9547f4e1 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C @@ -327,23 +327,25 @@ bool Foam::twoWayMPI::couple() const } // models with exact timing exists + label commandLines(0); if(exactTiming) { // extension for more liggghtsCommands active the same time: // sort interrupt list within this run period // keep track of corresponding liggghtsCommand int DEMstepsRun(0); + forAll(interruptTimes,j) { // set run command till interrupt DEMstepsRun += DEMstepsToInterrupt[j]; particleCloud_.liggghtsCommand()[runComNr]().set(DEMstepsToInterrupt[j]); - const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(); + const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0); Info << "Executing run command: '"<< command <<"'"<< endl; lmp->input->one(command); // run liggghts command with exact timing - command = particleCloud_.liggghtsCommand()[lcModel[j]]().command(); + command = particleCloud_.liggghtsCommand()[lcModel[j]]().command(0); Info << "Executing command: '"<< command <<"'"<< endl; lmp->input->one(command); } @@ -352,7 +354,7 @@ bool Foam::twoWayMPI::couple() const if(particleCloud_.liggghtsCommand()[runComNr]().runCommand(couplingStep())) { particleCloud_.liggghtsCommand()[runComNr]().set(couplingInterval() - DEMstepsRun); - const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(); + const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0); Info << "Executing run command: '"<< command <<"'"<< endl; lmp->input->one(command); } @@ -366,9 +368,13 @@ bool Foam::twoWayMPI::couple() const particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep()) ) { - const char* command = particleCloud_.liggghtsCommand()[i]().command(); - Info << "Executing command: '"<< command <<"'"<< endl; - lmp->input->one(command); + commandLines=particleCloud_.liggghtsCommand()[i]().commandLines(); + for(int j=0;jinput->one(command); + } } } } @@ -379,9 +385,13 @@ bool Foam::twoWayMPI::couple() const { if(particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep())) { - const char* command = particleCloud_.liggghtsCommand()[i]().command(); - Info << "Executing command: '"<< command <<"'"<< endl; - lmp->input->one(command); + commandLines=particleCloud_.liggghtsCommand()[i]().commandLines(); + for(int j=0;jinput->one(command); + } } } } diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H index 0afe3977..2efc6054 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H @@ -167,7 +167,7 @@ public: word myType() const{return typeName; }; - void setCG() const { particleCloud_.forceM(0).setCG(lmp->force->cg()); }; + void setCG() const { particleCloud_.setCG(lmp->force->cg()); }; }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C index abf2ef58..a8ef778a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C @@ -71,7 +71,12 @@ DiFeliceDrag::DiFeliceDrag rho_(sm.mesh().lookupObject (densityFieldName_)), voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), - interpolation_(false) + interpolation_(false), + splitImplicitExplicit_(false), + UsFieldName_(propsDict_.lookup("granVelFieldName")), + UsField_(sm.mesh().lookupObject (UsFieldName_)), + scaleDia_(1.), + scaleDrag_(1.) { //Append the field names to be probed particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat"); @@ -89,7 +94,19 @@ DiFeliceDrag::DiFeliceDrag Info << "using interpolated value of U." << endl; interpolation_=true; } - particleCloud_.checkCG(false); + if (propsDict_.found("splitImplicitExplicit")) + { + Info << "will split implicit / explicit force contributions." << endl; + splitImplicitExplicit_ = true; + if(!interpolation_) + Info << "WARNING: will only consider fluctuating particle velocity in implicit / explicit force split!" << endl; + } + + particleCloud_.checkCG(true); + if (propsDict_.found("scale")) + scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + if (propsDict_.found("scaleDrag")) + scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); } @@ -103,6 +120,13 @@ DiFeliceDrag::~DiFeliceDrag() void DiFeliceDrag::setForce() const { + if (scaleDia_ > 1) + Info << "DiFeliceDrag using scale = " << scaleDia_ << endl; + else if (particleCloud_.cg() > 1){ + scaleDia_=particleCloud_.cg(); + Info << "DiFeliceDrag using scale from liggghts cg = " << scaleDia_ << endl; + } + // get viscosity field #ifdef comp const volScalarField nufField = particleCloud_.turbulence().mu() / rho_; @@ -123,6 +147,11 @@ void DiFeliceDrag::setForce() const scalar magUr(0); scalar Rep(0); scalar Cd(0); + + vector UfluidFluct(0,0,0); + vector UsFluct(0,0,0); + vector dragExplicit(0,0,0); + scalar dragCoefficient(0); interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); @@ -136,7 +165,7 @@ void DiFeliceDrag::setForce() const cellI = particleCloud_.cellIDs()[index][0]; drag = vector(0,0,0); - + if (cellI > -1) // particle Found { if(interpolation_) @@ -158,11 +187,13 @@ void DiFeliceDrag::setForce() const magUr = mag(Ur); Rep = 0; Cd = 0; + dragCoefficient = 0; if (magUr > 0) { + // calc particle Re Nr - Rep = ds*voidfraction*magUr/(nuf+SMALL); + Rep = ds/scaleDia_*voidfraction*magUr/(nuf+SMALL); // calc fluid drag Coeff Cd = sqr(0.63 + 4.8/sqrt(Rep)); @@ -171,24 +202,45 @@ void DiFeliceDrag::setForce() const scalar Xi = 3.7 - 0.65 * exp(-sqr(1.5-log10(Rep))/2); // calc particle's drag - drag = 0.125*Cd*rho*M_PI*ds*ds*pow(voidfraction,(2-Xi))*magUr*Ur; - + dragCoefficient = 0.125*Cd*rho + *M_PI + *ds*ds + *scaleDia_ + *pow(voidfraction,(2-Xi))*magUr + *scaleDrag_; if (modelType_=="B") - drag /= voidfraction; + dragCoefficient /= voidfraction; + + drag = dragCoefficient*Ur; //total drag force! + + //Split forces + if(splitImplicitExplicit_) + { + UfluidFluct = Ufluid - U_[cellI]; + UsFluct = Us - UsField_[cellI]; + dragExplicit = dragCoefficient*(UfluidFluct - UsFluct); //explicit part of force + } + } - if(verbose_ && index >100 && index <102) + if(verbose_ && index >-1 && index <102) { Pout << "index = " << index << endl; Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; - Pout << "ds = " << ds << endl; + Pout << "ds/scale = " << ds/scaleDia_ << endl; Pout << "rho = " << rho << endl; Pout << "nuf = " << nuf << endl; Pout << "voidfraction = " << voidfraction << endl; Pout << "Rep = " << Rep << endl; Pout << "Cd = " << Cd << endl; - Pout << "drag = " << drag << endl; + Pout << "drag (total) = " << drag << endl; + if(splitImplicitExplicit_) + { + Pout << "UfluidFluct = " << UfluidFluct << endl; + Pout << "UsFluct = " << UsFluct << endl; + Pout << "dragExplicit = " << dragExplicit << endl; + } } //Set value fields and write the probe @@ -205,7 +257,15 @@ void DiFeliceDrag::setForce() const } // set force on particle if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] += drag[j]; - else for(int j=0;j<3;j++) impForces()[index][j] += drag[j]; + else //implicit treatment, taking explicit force contribution into account + { + for(int j=0;j<3;j++) + { + impForces()[index][j] += drag[j] - dragExplicit[j]; //only consider implicit part! + expForces()[index][j] += dragExplicit[j]; + } + } + for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j]; } //} diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H index 45a285a6..5d788bf6 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H @@ -76,6 +76,16 @@ private: bool interpolation_; // use interpolated U field values + bool splitImplicitExplicit_; // use splitting of implicit and explict force contribution + + word UsFieldName_; + + const volVectorField& UsField_; // the average particle velocity field (for implicit/expliti force split) + + mutable scalar scaleDia_; + + mutable scalar scaleDrag_; + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C index 0e6f34c0..c4f8e888 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C @@ -71,21 +71,34 @@ GidaspowDrag::GidaspowDrag voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), phi_(readScalar(propsDict_.lookup("phi"))), - interpolation_(false) + interpolation_(false), + scaleDia_(1.), + scaleDrag_(1.) { //Append the field names to be probed particleCloud_.probeM().initialize(typeName, "gidaspowDrag.logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().scalarFields_.append("Rep"); - particleCloud_.probeM().scalarFields_.append("beta"); + particleCloud_.probeM().scalarFields_.append("betaP"); particleCloud_.probeM().scalarFields_.append("voidfraction"); particleCloud_.probeM().writeHeader(); if (propsDict_.found("verbose")) verbose_=true; if (propsDict_.found("treatExplicit")) treatExplicit_=true; if (propsDict_.found("interpolation")) interpolation_=true; - particleCloud_.checkCG(false); + if (propsDict_.found("implDEM")) + { + treatExplicit_=false; + implDEM_=true; + setImpDEMdrag(); + Info << "Using implicit DEM drag formulation." << endl; + } + particleCloud_.checkCG(true); + if (propsDict_.found("scale")) + scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + if (propsDict_.found("scaleDrag")) + scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); Info << "Gidaspow - interpolation switch: " << interpolation_ << endl; } @@ -101,6 +114,13 @@ GidaspowDrag::~GidaspowDrag() void GidaspowDrag::setForce() const { + if (scaleDia_ > 1) + Info << "Gidaspow using scale = " << scaleDia_ << endl; + else if (particleCloud_.cg() > 1){ + scaleDia_=particleCloud_.cg(); + Info << "Gidaspow using scale from liggghts cg = " << scaleDia_ << endl; + } + // get viscosity field #ifdef comp const volScalarField nufField = particleCloud_.turbulence().mu() / rho_; @@ -111,6 +131,8 @@ void GidaspowDrag::setForce() const vector position(0,0,0); scalar voidfraction(1); vector Ufluid(0,0,0); + vector drag(0,0,0); + label cellI=0; vector Us(0,0,0); vector Ur(0,0,0); @@ -122,9 +144,9 @@ void GidaspowDrag::setForce() const scalar Vs(0); scalar localPhiP(0); - scalar CdMagUrLag(0); //Cd of the very particle - scalar KslLag(0); //momentum exchange of the very particle (per unit volume) - scalar beta(0); //momentum exchange of the very particle + scalar CdMagUrLag(0); //Cd of the very particle + scalar KslLag(0); //momentum exchange of the very particle (per unit volume) + scalar betaP(0); //momentum exchange of the very particle interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); @@ -135,8 +157,12 @@ void GidaspowDrag::setForce() const { //if(mask[index][0]) //{ - vector drag(0,0,0); - label cellI = particleCloud_.cellIDs()[index][0]; + cellI = particleCloud_.cellIDs()[index][0]; + drag = vector(0,0,0); + betaP = 0; + Vs = 0; + Ufluid =vector(0,0,0); + voidfraction=0; if (cellI > -1) // particle Found { @@ -148,8 +174,8 @@ void GidaspowDrag::setForce() const Ufluid = UInterpolator_.interpolate(position,cellI); //Ensure interpolated void fraction to be meaningful // Info << " --> voidfraction: " << voidfraction << endl; - if(voidfraction>1.00) voidfraction = 1.0f; - if(voidfraction<0.10) voidfraction = 0.10f; + if(voidfraction>1.00) voidfraction = 1.0; + if(voidfraction<0.10) voidfraction = 0.10; } else { @@ -171,30 +197,30 @@ void GidaspowDrag::setForce() const //Compute specific drag coefficient (i.e., Force per unit slip velocity and per m³ SUSPENSION) if(voidfraction > 0.8) //dilute { - Rep=ds*voidfraction*magUr/nuf; - CdMagUrLag = (24.0*nuf/(ds*voidfraction)) //1/magUr missing here, but compensated in expression for KslLag! + Rep=ds/scaleDia_*voidfraction*magUr/nuf; + CdMagUrLag = (24.0*nuf/(ds/scaleDia_*voidfraction)) //1/magUr missing here, but compensated in expression for KslLag! *(scalar(1)+0.15*Foam::pow(Rep, 0.687)); KslLag = 0.75*( rho*localPhiP*voidfraction*CdMagUrLag / - (ds*Foam::pow(voidfraction,2.65)) + (ds/scaleDia_*Foam::pow(voidfraction,2.65)) ); } else //dense { KslLag = (150*Foam::pow(localPhiP,2)*nuf*rho)/ - (voidfraction*ds*ds+SMALL) + (voidfraction*ds/scaleDia_*ds/scaleDia_+SMALL) + (1.75*(localPhiP) * magUr * rho)/ - ((ds)); + ((ds/scaleDia_)); } // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE) - beta = KslLag / localPhiP; + betaP = KslLag / localPhiP; // calc particle's drag - drag = Vs * beta * Ur; + drag = Vs * betaP * Ur * scaleDrag_; if (modelType_=="B") drag /= voidfraction; @@ -206,11 +232,13 @@ void GidaspowDrag::setForce() const Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; Pout << "ds = " << ds << endl; + Pout << "ds/scale = " << ds/scaleDia_ << endl; + Pout << "phi = " << phi_ << endl; Pout << "rho = " << rho << endl; Pout << "nuf = " << nuf << endl; Pout << "voidfraction = " << voidfraction << endl; Pout << "Rep = " << Rep << endl; - Pout << "beta = " << beta << endl; + Pout << "betaP = " << betaP << endl; Pout << "drag = " << drag << endl; } @@ -221,7 +249,7 @@ void GidaspowDrag::setForce() const vValues.append(drag); //first entry must the be the force vValues.append(Ur); sValues.append(Rep); - sValues.append(beta); + sValues.append(betaP); sValues.append(voidfraction); particleCloud_.probeM().writeProbe(index, sValues, vValues); } @@ -230,7 +258,21 @@ void GidaspowDrag::setForce() const // set force on particle if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] += drag[j]; else for(int j=0;j<3;j++) impForces()[index][j] += drag[j]; - for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j]; + + // set Cd + if(implDEM_) + { + for(int j=0;j<3;j++) fluidVel()[index][j]=Ufluid[j]; + + if (modelType_=="B" && cellI > -1) + Cds()[index][0] = Vs*betaP/voidfraction*scaleDrag_; + else + Cds()[index][0] = Vs*betaP*scaleDrag_; + + }else{ + for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j]; + } + //}// end if mask }// end loop particles } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H index f9c881b9..978ca691 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H @@ -80,6 +80,10 @@ private: bool interpolation_; // use interpolated field values + mutable scalar scaleDia_; + + mutable scalar scaleDrag_; + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C index 50f36bef..a6021394 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C @@ -115,8 +115,8 @@ void KochHillDrag::setForce() const { if (scaleDia_ > 1) Info << "KochHill using scale = " << scaleDia_ << endl; - else if (cg() > 1){ - scaleDia_=cg(); + else if (particleCloud_.cg() > 1){ + scaleDia_=particleCloud_.cg(); Info << "KochHill using scale from liggghts cg = " << scaleDia_ << endl; } @@ -158,6 +158,7 @@ void KochHillDrag::setForce() const betaP = 0; Vs = 0; Ufluid =vector(0,0,0); + voidfraction=0; if (cellI > -1) // particle Found { @@ -233,6 +234,7 @@ void KochHillDrag::setForce() const Pout << "nuf = " << nuf << endl; Pout << "voidfraction = " << voidfraction << endl; Pout << "Rep = " << Rep << endl; + Pout << "betaP = " << betaP << endl; Pout << "drag = " << drag << endl; } @@ -257,10 +259,10 @@ void KochHillDrag::setForce() const { for(int j=0;j<3;j++) fluidVel()[index][j]=Ufluid[j]; - if (modelType_=="B") - Cds()[index][0] = Vs*betaP/voidfraction; + if (modelType_=="B" && cellI > -1) + Cds()[index][0] = Vs*betaP/voidfraction*scaleDrag_; else - Cds()[index][0] = Vs*betaP; + Cds()[index][0] = Vs*betaP*scaleDrag_; }else{ for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j]; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C index 0ad0b471..f68ecba0 100755 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C @@ -140,12 +140,10 @@ void KochHillRWDrag::setForce() const // realloc the arrays reAllocArrays(); -Info << "huhu" << endl; - if (scale_ > 1) Info << "KochHillRW using scale = " << scale_ << endl; - else if (cg() > 1){ - scale_=cg(); + else if (particleCloud_.cg() > 1){ + scale_=particleCloud_.cg(); Info << "KochHillRW using scale from liggghts cg = " << scale_ << endl; } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/checkCouplingInterval/checkCouplingInterval.C b/src/lagrangian/cfdemParticle/subModels/forceModel/checkCouplingInterval/checkCouplingInterval.C index 6d61e0e2..97ad879f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/checkCouplingInterval/checkCouplingInterval.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/checkCouplingInterval/checkCouplingInterval.C @@ -99,7 +99,7 @@ void checkCouplingInterval::setForce() const if (cellI > -1) // particle Found { - scaledRad = particleCloud_.radius(index)/cg(); + scaledRad = particleCloud_.radius(index)/particleCloud_.cg(); tauP = rhoP_*4*scaledRad*scaledRad/ (18 * nufField[cellI] * rho_[cellI]); minTauP = min(minTauP,tauP); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.C b/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.C index c8b75e20..3e14866f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.C @@ -74,9 +74,8 @@ fieldTimeAverage::fieldTimeAverage // create time average scalar fields scalarFields_.setSize(scalarFieldNames_.size()); - if (propsDict_.found("startTime")){ + if (propsDict_.found("startTime")) startTime_=readScalar(propsDict_.lookup("startTime")); - } for (int i=0;i < scalarFieldNames_.size(); i++) { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.H b/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.H index 5de898ec..f9e05696 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/fieldTimeAverage/fieldTimeAverage.H @@ -73,8 +73,7 @@ private: mutable PtrList vectorFields_; mutable double nrAverages_; - - + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C index d4b94e37..d87a304a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C @@ -85,8 +85,8 @@ forceModel::forceModel ), coupleForce_(true), modelType_(sm.modelType()), - cg_(1.), - probeIt_(sm.probeM().active()) + probeIt_(sm.probeM().active()), + requiresEx_(false) {} @@ -130,6 +130,21 @@ void forceModel::manipulateScalarField(volScalarField& field) const Info << "no scalar manipulation done" << endl; // do nothing } + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void forceModel::repartitionImExForces() const +{ + if(particleCloud_.imExSplitFactor()<1.0) + { + Info << "Will re-partition split of implicit and explicit forces: alpha = " + << particleCloud_.imExSplitFactor() << endl; + // Update implicit particle + expParticleForces_ += (1.0-particleCloud_.imExSplitFactor())*impParticleForces_; + impParticleForces_ *= particleCloud_.imExSplitFactor(); + } +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H index 82ebb0ad..8a474797 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H @@ -75,9 +75,9 @@ protected: const word modelType_; - mutable scalar cg_; - bool probeIt_; + + bool requiresEx_; //requires a orientation vector public: @@ -150,17 +150,12 @@ public: inline const bool& coupleForce() const { return coupleForce_;}; - inline const scalar cg() const { return cg_; }; - - void const setCG(double cg) const - { - cg_ = cg; - Info << "cg is set to: " << cg_ << endl; - }; - void const setImpDEMdrag() const {particleCloud_.impDEMdrag_=true;}; - - + + virtual inline bool& requiresEx() { return requiresEx_;}; + + //Repartition Implixit/Explicit forces + void repartitionImExForces() const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C index edd42807..e3073a21 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C @@ -100,7 +100,15 @@ void noDrag::setForce() const if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] = 0.; else for(int j=0;j<3;j++) impForces()[index][j] = 0.; } - if(noDEMForce_)for(int j=0;j<3;j++) DEMForces()[index][j] = 0.; + if(noDEMForce_) + { + for(int j=0;j<3;j++) DEMForces()[index][j] = 0.; + if(particleCloud_.impDEMdrag()) + { + Cds()[index][0] = 0.; + for(int j=0;j<3;j++) fluidVel()[index][j] = 0.; + } + } } } } diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C index 99ac99c2..c13c59e6 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C @@ -89,73 +89,9 @@ execute::execute // check if verbose if (propsDict_.found("verbose")) verbose_=true; - bool addBlank = true; // std no blanks after each word - fileName add; - label numberCount=0; // nr of scalars inserted to command - label labelCount=0; // nr of labels inserted to command - - forAll(commandList_,i) - { - add = word(commandList_[i]); - - //- handle symbols - if (add == "$couplingInterval") - { - char h[50]; - sprintf(h,"%d",particleCloud_.dataExchangeM().couplingInterval()); - add = h; - } - else if (add=="dot") add = "."; - else if (add=="dotdot") add = ".."; - else if (add=="slash") add = "/"; - else if (add=="noBlanks") // no blanks after the following words - { - add = ""; - addBlank = false; - }else if (add=="blanks") // add a blank here and after the following words - { - add = ""; - addBlank = true; - }else if (add=="timeStamp") // next command will be a number read from scalarList_ - { - add = ""; - timeStamp_=true; - }else if (add=="number") // next command will be a number read from scalarList_ - { - if (!propsDict_.found("scalars")) - { - FatalError<<"you want to use a number in the command\n - specify a scalar list with all numbers" - << abort(FatalError); - - } - char h[50]; - sprintf(h,"%f",scalarList_[numberCount]); - add = h; - numberCount ++; - }else if (add=="label") // next command will be a number read from scalarList_ - { - if (!propsDict_.found("labels")) - { - FatalError<<"you want to use a label in the command\n - specify a label list with all numbers" - << abort(FatalError); - - } - char h[50]; - sprintf(h,"%d",labelList_[labelCount]); - add = h; - labelCount ++; - } - - // compose command - if (addBlank) - { - command_ += add + " "; - }else - { - command_ += add; - } - } + parseCommandList(commandList_, labelList_, scalarList_, command_, propsDict_, timeStamp_); Info << "liggghtsCommand " << command_ << endl; + strCommand_=string(command_); checkTimeMode(propsDict_); @@ -172,7 +108,7 @@ execute::~execute() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const char* execute::command() +const char* execute::command(int commandLine) { return strCommand_.c_str(); } diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.H b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.H index bc9111e8..1b73e65a 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.H +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.H @@ -98,7 +98,7 @@ public: word name(){return myName_;}; - const char* command(); + const char* command(int); bool runCommand(int); diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C index a5d78a42..ae85ed1c 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C @@ -71,6 +71,7 @@ liggghtsCommandModel::liggghtsCommandModel lastCouplingStep_(-1), couplingStepInterval_(0), exactTiming_(false), + commandLines_(1), verbose_(false) {} @@ -236,6 +237,74 @@ DynamicList liggghtsCommandModel::executionsWithinPeriod(scalar TSstart, return executions; } + +void liggghtsCommandModel::parseCommandList(wordList& commandList,labelList& labelList,scalarList& scalarList,word& command, dictionary& propsDict, bool timeStamp) +{ + bool addBlank = true; // std no blanks after each word + fileName add; + label numberCount=0; // nr of scalars inserted to command + label labelCount=0; // nr of labels inserted to command + + forAll(commandList,i) + { + add = word(commandList[i]); + + //- handle symbols + if (add == "$couplingInterval") + { + char h[50]; + sprintf(h,"%d",particleCloud_.dataExchangeM().couplingInterval()); + add = h; + } + else if (add=="dot") add = "."; + else if (add=="dotdot") add = ".."; + else if (add=="slash") add = "/"; + else if (add=="noBlanks") // no blanks after the following words + { + add = ""; + addBlank = false; + }else if (add=="blanks") // add a blank here and after the following words + { + add = ""; + addBlank = true; + }else if (add=="timeStamp") // next command will be a number read from labelList + { + add = ""; + timeStamp=true; + }else if (add=="number") // next command will be a number read from labelList + { + /*if (!propsDict.found("scalars")) + { + FatalError<<"you want to use a number in the command\n - specify a scalar list with all numbers" + << abort(FatalError); + }*/ + char h[50]; + sprintf(h,"%f",scalarList[numberCount]); + add = h; + numberCount ++; + }else if (add=="label") // next command will be a number read from labelList + { + /*if (!propsDict.found("labels")) + { + FatalError<<"you want to use a label in the command\n - specify a label list with all numbers" + << abort(FatalError); + }*/ + char h[50]; + sprintf(h,"%d",labelList[labelCount]); + add = h; + labelCount ++; + } + + // compose command + if (addBlank) + { + command += add + " "; + }else + { + command += add; + } + } +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.H b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.H index 18ffcdd4..cba32881 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.H +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.H @@ -89,6 +89,8 @@ protected: bool exactTiming_; + label commandLines_; + bool verbose_; public: @@ -142,7 +144,7 @@ public: // Member Functions - virtual const char* command()=0; + virtual const char* command(int)=0; void checkTimeMode(dictionary&); @@ -168,6 +170,9 @@ public: bool exactTiming(){return exactTiming_;}; + label commandLines(){return commandLines_;}; + + void parseCommandList(wordList&, labelList&, scalarList&, word&, dictionary&, bool=false); }; diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.C index bf7d62d5..3605ce97 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.C @@ -112,7 +112,7 @@ readLiggghtsData::~readLiggghtsData() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const char* readLiggghtsData::command() +const char* readLiggghtsData::command(int commandLine) { char h[50]; sprintf(h,"_%d",insertionNr_); diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.H b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.H index fd257244..345d2c4d 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.H +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.H @@ -93,7 +93,7 @@ public: word name(){return myName_;}; - const char* command(); + const char* command(int); bool runCommand(int); diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.C index 17a23b3a..2969a052 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.C @@ -93,7 +93,7 @@ runLiggghts::~runLiggghts() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const char* runLiggghts::command() +const char* runLiggghts::command(int commandLine) { return strCommand_.c_str(); } diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.H b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.H index e226c84b..39aaa3ee 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.H +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.H @@ -88,7 +88,7 @@ public: // Member Functions - const char* command(); + const char* command(int); string createCommand(word="",int=0,word="",word="",word="",word=""); diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C index 91214a2f..4d639585 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C @@ -115,7 +115,7 @@ writeLiggghts::~writeLiggghts() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const char* writeLiggghts::command() +const char* writeLiggghts::command(int commandLine) { return strCommand_.c_str(); } diff --git a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.H b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.H index d1bbf3bd..c6e1fdc6 100644 --- a/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.H +++ b/src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.H @@ -93,7 +93,7 @@ public: // Member Functions - const char* command(); + const char* command(int); bool runCommand(int); diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U index bab71f9e..d6975143 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U @@ -51,29 +51,12 @@ boundaryField value uniform (0 0 0);*/ // 2.1.x + //type uniformFixedValue; type uniformFixedValueVoidfraction; uniformValue table ( (0.000 (0 0 0.002)) - (0.020 (0 0 0.002)) - (0.021 (0 0 0.004)) - (0.040 (0 0 0.004)) - (0.041 (0 0 0.006)) - (0.060 (0 0 0.006)) - (0.061 (0 0 0.008)) - (0.080 (0 0 0.008)) - (0.081 (0 0 0.010)) - (0.100 (0 0 0.010)) - (0.101 (0 0 0.012)) - (0.120 (0 0 0.012)) - (0.121 (0 0 0.014)) - (0.140 (0 0 0.014)) - (0.141 (0 0 0.016)) - (0.160 (0 0 0.016)) - (0.161 (0 0 0.018)) - (0.180 (0 0 0.018)) - (0.181 (0 0 0.020)) - (0.200 (0 0 0.020)) + (0.100 (0 0 0.030)) ); } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p index d111e6be..be1cbf5c 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p @@ -35,8 +35,9 @@ boundaryField outlet { //type zeroGradient; + type fixedValue; - value uniform 100000; + value uniform 1.0e5; } } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho index 2c9e8cd0..3dc9d41a 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho @@ -16,7 +16,7 @@ FoamFile dimensions [1 -3 0 0 0 0 0]; -internalField uniform 30; +internalField uniform 10; boundaryField { diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction index 8298f8c4..ac2147e2 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction @@ -28,10 +28,14 @@ boundaryField outlet { type zeroGradient; + //type fixedValue; + //value uniform 1; } inlet { type zeroGradient; + //type fixedValue; + //value uniform 1; } } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties index 3bc089ba..5764aac3 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties @@ -25,42 +25,44 @@ FoamFile //===========================================================================// // sub-models & settings +solveFlow true; modelType "A"; // A or B couplingInterval 100; -voidFractionModel divided; +voidFractionModel divided;//centre;// -locateModel engine; +locateModel engine;//turboEngineM2M;// meshMotionModel noMeshMotion; regionModel allRegion; -IOModel "basicIO"; +IOModel basicIO; probeModel off; -dataExchangeModel twoWayMPI;//twoWayFiles;//oneWayVTK;// +dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;// averagingModel dense;//dilute;// clockModel off;//standardClock;// -smoothingModel off; +smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; // forceModels ( - //GidaspowDrag - //DiFeliceDrag - KochHillDrag + //KochHillDrag //BeetstraDrag + //noDrag gradPForce viscForce + GidaspowDrag //Archimedes //volWeightedAverage //totalMomentumExchange + particleCellVolume ); momCoupleModels @@ -73,6 +75,21 @@ turbulenceModelType "RASProperties";//"LESProperties";// //===========================================================================// // sub-model properties +localPSizeDiffSmoothingProps +{ + lowerLimit 0.1; + upperLimit 1e10; + dSmoothingLength 1.5e-3; + Csmoothing 1.0; +} + +constDiffSmoothingProps +{ + lowerLimit 0.1; + upperLimit 1e10; + smoothingLength 1.5e-3; +} + implicitCoupleProps { velFieldName "U"; @@ -91,7 +108,7 @@ gradPForceProps densityFieldName "rho"; voidfractionFieldName "voidfraction"; velocityFieldName "U"; - //interpolation; + interpolation; } viscForceProps @@ -100,6 +117,11 @@ viscForceProps densityFieldName "rho"; interpolation; } +noDragProps +{ + //noDEMForce; + //keepCFDForce; +} volWeightedAverageProps { scalarFieldNames @@ -123,22 +145,24 @@ totalMomentumExchangeProps } GidaspowDragProps { - velFieldName "U"; - densityFieldName "rho"; -} -DiFeliceDragProps -{ + verbose; velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; + phi 1; + //interpolation; // this case does not like interpolation + implDEM_; } KochHillDragProps { + //verbose; velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; - //verbose; + //interpolation; + implDEM_; + verbose; } BeetstraDragProps @@ -162,6 +186,13 @@ virtualMassForceProps densityFieldName "rho"; } +particleCellVolumeProps +{ + upperThreshold 0.999; + lowerThreshold 0.; + verbose; +} + oneWayVTKProps { couplingFilename "vtk_out%4.4d.vtk"; @@ -175,7 +206,7 @@ twoWayFilesProps centreProps { - alphaMin 0.10; + alphaMin 0.1; } engineProps @@ -185,7 +216,7 @@ engineProps dividedProps { - alphaMin 0.1; + alphaMin 0.05; scaleUpVol 1.0; } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m index 9b2c22ff..1532e2ef 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/octave/totalPressureDrop.m @@ -5,7 +5,7 @@ clc; %====================================% % simulation data 1 %====================================% -rhoG = 30 % density in kg/m3 +rhoG = 10 % density in kg/m3 %path = '../probes/0/p'; path = '../postProcessing/probes/0/p'; columns=22; @@ -13,7 +13,7 @@ headerlines=4; data = loaddata(path,columns,headerlines); data=transpose(data); [x,y]=size(data) -dp_sim = (data(:,2)-data(:,y))*rhoG; +dp_sim = (data(:,2)-data(:,y))*rhoG; %conversion to Pa! t_sim = data(:,1); %fprintf('final pressureDrop of sim = %f Pa\n',dp_sim(length(dp_sim)) ) @@ -25,18 +25,18 @@ t_sim = data(:,1); % Ergun Equation %=================== fprintf('\ncalc Ergun eqn:\n') -dp = 0.0008 % particle diameter +dp = 0.0005 % particle diameter phip = 1 % sphericity -epsilon = 0.451335 % void fraction +epsilon = 1-0.518497 % void fraction NOTE: varies with voidfraction!!! Ustart = 0.002 -Uend = 0.02 +Uend = 0.03 timeStepSize = 0.001; % time interval of pressure data Tstart = 0; Tend = t_sim(length(t_sim)); deltaU=(Uend-Ustart)/((Tend-Tstart)/timeStepSize); U = Ustart+deltaU:deltaU:Uend; % velocity over time Ua = U / epsilon; % physical velocity -L = 0.0156 % length of bed +L = 0.0135 % length of bed nuG = 1.5*10^-4 % kinemat Visk in m2/s muG = nuG*rhoG % dynam visc in Pa s @@ -55,14 +55,31 @@ fprintf('so the result does not depend on density\n') %================================== rhoP = 7000 % particle density in kg/m3 g = 9.81 % gravity m/s2 + +% WARNING!!! +% these equations can give wrong results if applied to the wrong regime!!! + Umf = dp^2*(rhoP-rhoG)*g/(150*muG)*(epsilon^3*phip^2)/(1-epsilon) -ReMF = Umf*dp*rhoG/muG % must be <20 !!! -%Umf = sqrt(phip*dp^2/1.75*(rhoP-rhoG)/rhoG*g*epsilon^3) % Re>1000 +ReMF = Umf*dp*rhoG/muG +#Umf = sqrt(dp*(rhoP-rhoG)*g/(1.75*rhoG)*epsilon^3*phip) +#ReMF = Umf*dp*rhoG/muG + + +if(ReMF<20) + fprintf('applying eqn1 for Umf.\n') +elseif(ReMF>20 && ReMF<1000) + fprintf('applying eqn1 for Umf.\n') +elseif (ReMF>=1000) + fprintf('applying eqn2 for Umf.\n') + Umf = sqrt(dp*(rhoP-rhoG)*g/(1.75*rhoG)*epsilon^3*phip); + ReMF = Umf*dp*rhoG/muG; +end dpUmf= L * ( 150*((1-epsilon)^2/epsilon^3)*((muG.*Umf)/(phip*dp)^2) +1.75*((1-epsilon)/epsilon^3)*((rhoG.*Umf.^2)/(phip*dp)) - ); + ) +dpUmf=(L*(1-epsilon)*(rhoP-rhoG)*g) %====================================% % plot data @@ -80,7 +97,7 @@ figure(1) plot(U,dpErgun,U,dp_sim,[Umf,Uend],dpUmf*ones(1,2)) title("Ergun pressure drop vs. simulation") a=strcat("analytical (Ergun), Umf=",num2str(Umf),", dpUmf=",num2str(dpUmf)); -legend(a,"simulation") +legend(a,"simulation","analyt. deltaP at Umf") xlabel("velocity in [m/s]") ylabel("pressure drop [Pa]") axis([0,Uend,0,dpErgun(length(dpErgun))]) diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict index 9e8fbc16..85e89dbe 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict @@ -23,13 +23,13 @@ startTime 0; stopAt endTime; -endTime 0.2;//0.01; +endTime 0.1; deltaT 0.001; writeControl adjustableRunTime; -writeInterval 0.005;//0.01; +writeInterval 0.01; purgeWrite 0; @@ -54,6 +54,7 @@ libs ( "libfiniteVolumeCFDEM.so" ); functions ( + probes { type probes; @@ -87,7 +88,7 @@ functions ); // Fields to be probed - fields ( p U voidfraction volAverage_voidfraction); + fields ( p U voidfraction volAverage_voidfraction voidfractionNext voidfractionPrev); // Write at same frequency as fields outputControl timeStep;//outputTime; diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict index bccec657..4de14078 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict @@ -30,7 +30,7 @@ numberOfSubdomains 4; simpleCoeffs { - n (1 2 2); + n (2 2 1); delta 0.001; } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution index 955f5b14..46c6e11d 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution @@ -17,7 +17,7 @@ FoamFile solvers { - p + "(p)" { solver PCG; preconditioner DIC; @@ -33,7 +33,7 @@ solvers relTol 0; } - U + "(U|k|epsilon|R|nuTilda)" { solver PBiCG; preconditioner DILU; @@ -41,34 +41,10 @@ solvers relTol 0; } - k + "(voidfraction|Us|Ksl|dSmoothing)" { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - epsilon - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - R - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - nuTilda - { - solver PBiCG; - preconditioner DILU; + solver PCG; + preconditioner DIC; tolerance 1e-05; relTol 0; } diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init index acb6874b..3fdff28d 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init @@ -7,9 +7,9 @@ echo both ####################################################### # variables # ####################################################### -variable cg equal 2 -variable r0 equal 0.0004 -variable nPorg equal 19531 +variable cg equal 1 +variable r0 equal 0.00025 +variable nPorg equal 71000 variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg}) ####################################################### @@ -22,6 +22,7 @@ boundary m m m newton off units si +processors 2 2 1 region reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box create_box 1 reg @@ -36,24 +37,23 @@ fix m1 all property/global youngsModulus peratomtype 5.e6 fix m2 all property/global poissonsRatio peratomtype 0.45 fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3 fix m4 all property/global coefficientFriction peratomtypepair 1 0.5 -fix m5 all property/global characteristicVelocity scalar 1. #pair style -pair_style gran model hooke tangential history +pair_style gran model hertz tangential history #Hertzian without cohesion pair_coeff * * #timestep, gravity -timestep 0.000005 +timestep 0.00001 fix gravi all gravity 9.81 vector 0.0 0.0 -1.0 #walls -fix zwalls1 all wall/gran model hooke tangential history primitive type 1 zplane 0.0 -fix zwalls2 all wall/gran model hooke tangential history primitive type 1 zplane 0.0553 -fix cylwalls all wall/gran model hooke tangential history primitive type 1 zcylinder 0.01385 0. 0. +fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0 +fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 0.0553 +fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.01385 0. 0. #particle distributions and insertion region bc cylinder z 0.0 0.0 0.012 0. 0.055 units box -fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 200 radius constant ${r0} +fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 7000 radius constant ${r0} fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0 fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -1. insert_every once overlapcheck yes all_in yes particles_in_region ${nPscal} region bc diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume index dd167a16..ea12ebe0 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume @@ -7,9 +7,9 @@ echo both ####################################################### # variables # ####################################################### -variable cg equal 2 -variable r0 equal 0.0004 -variable nPorg equal 19531 +variable cg equal 3 +variable r0 equal 0.00025 +variable nPorg equal 71000 variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg}) ####################################################### @@ -22,28 +22,24 @@ boundary m m m newton off units si -processors 1 2 2 +processors 2 2 1 #read the restart file read_restart ../DEM/liggghts.restart_coarseGraining_${cg} -#do not do this here, the simulation box is in the restart file! -#region reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box -#create_box 1 reg - -neighbor 0.0005 bin +neighbor 0.001 bin neigh_modify delay 0 + #Material properties required for new pair styles fix m1 all property/global youngsModulus peratomtype 5.e6 fix m2 all property/global poissonsRatio peratomtype 0.45 fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3 fix m4 all property/global coefficientFriction peratomtypepair 1 0.5 -fix m5 all property/global characteristicVelocity scalar 1. #pair style -pair_style gran model hooke tangential history #Hertzian without cohesion +pair_style gran model hertz tangential history #Hertzian without cohesion pair_coeff * * #timestep, gravity @@ -51,16 +47,14 @@ timestep 0.00001 fix gravi all gravity 9.81 vector 0.0 0.0 -1.0 #walls -fix zwalls1 all wall/gran model hooke tangential history primitive type 1 zplane 0.0 -fix zwalls2 all wall/gran model hooke tangential history primitive type 1 zplane 0.0553 -fix cylwalls all wall/gran model hooke tangential history primitive type 1 zcylinder 0.01385 0. 0. +fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0 +fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 0.0553 +fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.01385 0. 0. -# change the particles density -set group all density 7000 #cfd coupling fix cfd all couple/cfd couple_every 100 mpi -fix cfd2 all couple/cfd/force +fix cfd2 all couple/cfd/force #/implicit #apply nve integration to all particles that are inserted as single particles fix integr all nve/sphere @@ -69,11 +63,11 @@ fix integr all nve/sphere compute centerOfMass all com #compute total dragforce -compute dragtotal all reduce sum f_dragforce[1] f_dragforce[2] f_dragforce[3] +#compute dragtotal all reduce sum f_dragforce[1] f_dragforce[2] f_dragforce[3] #screen output compute 1 all erotate/sphere -thermo_style custom step atoms ke c_1 vol c_centerOfMass[3] c_dragtotal[1] c_dragtotal[2] c_dragtotal[3] +thermo_style custom step atoms ke c_1 vol c_centerOfMass[3] #c_dragtotal[1] c_dragtotal[2] c_dragtotal[3] thermo 10 thermo_modify lost ignore norm no compute_modify thermo_temp dynamic yes @@ -81,7 +75,7 @@ compute_modify thermo_temp dynamic yes #insert the first particles so that dump is not empty dump myDump all stl 1 post/dump_*.stl #run 1 -dump dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius +dump dmp all custom 1000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz radius #f_dragforce[1] f_dragforce[2] f_dragforce[3] undump myDump run 1 diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy index ebf63dfd..e69de29b 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy @@ -1 +0,0 @@ -dummyfile diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh b/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh index 3ae20ab9..75569488 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh @@ -21,7 +21,7 @@ logfileName="log_$headerText" solverName="cfdemSolverPiso" nrProcs="4" machineFileName="none" # yourMachinefileName | none -debugMode="off" # on | off | strict +debugMode="off" # on | off| strict testHarnessPath="$CFDEM_TEST_HARNESS_PATH" runOctave="true" postproc="false" @@ -82,21 +82,10 @@ fi #- clean up case echo "deleting data at: $casePath :\n" -rm -r $casePath/CFD/0.* -rm -r $casePath/CFD/log.* -rm -r $casePath/CFD/octave/octave-core -rm -r $casePath/CFD/VTK -rm -r $casePath/CFD/processor* -rm -r $casePath/CFD/couplingFiles/* -rm -r $casePath/DEM/post/* -rm -r $casePath/DEM/log.* -rm -r $casePath/DEM/liggghts.restartCFDEM* -rm -r $casePath/CFD/probes* -rm -r $casePath/CFD/postProcessing -rm -r $casePath/CFD/particles -rm -r $casePath/CFD/lagrangian +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanCase rm -r $casePath/CFD/clockData +rm -r $casePath/DEM/post/* +(cd $casePath/DEM/post && touch dummy) echo "done" - -#- preserve post directory -echo "dummyfile" >> $casePath/DEM/post/dummy diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U index 9bd6c981..5fc0e4e2 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U @@ -55,24 +55,6 @@ boundaryField uniformValue table ( (0.000 (0 0 0.002)) - /*(0.010 (0 0 0.002)) - (0.011 (0 0 0.004)) - (0.020 (0 0 0.004)) - (0.021 (0 0 0.006)) - (0.030 (0 0 0.006)) - (0.031 (0 0 0.008)) - (0.040 (0 0 0.008)) - (0.041 (0 0 0.010)) - (0.050 (0 0 0.010)) - (0.051 (0 0 0.012)) - (0.060 (0 0 0.012)) - (0.061 (0 0 0.014)) - (0.070 (0 0 0.014)) - (0.071 (0 0 0.016)) - (0.080 (0 0 0.016)) - (0.081 (0 0 0.018)) - (0.090 (0 0 0.018)) - (0.091 (0 0 0.020))*/ (0.100 (0 0 0.020)) ); } diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties index f5336a0a..aa9fa97a 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties @@ -62,6 +62,7 @@ forceModels //volWeightedAverage //totalMomentumExchange //particleCellVolume + //fieldTimeAverage ); momCoupleModels @@ -107,7 +108,7 @@ gradPForceProps densityFieldName "rho"; voidfractionFieldName "voidfraction"; velocityFieldName "U"; - //interpolation; + interpolation; } viscForceProps @@ -143,6 +144,7 @@ GidaspowDragProps velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; + interpolation; phi 1; } DiFeliceDragProps @@ -158,6 +160,7 @@ KochHillDragProps velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; + interpolation; } BeetstraDragProps @@ -201,7 +204,7 @@ twoWayFilesProps centreProps { - alphaMin 0.10; + alphaMin 0.1; } engineProps diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/liggghtsCommands index 4ecd5aa5..a5ffcaed 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/liggghtsCommands +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/liggghtsCommands @@ -25,7 +25,7 @@ FoamFile liggghtsCommandModels ( - runLiggghts + runLiggghts ); // ************************************************************************* // diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m index 233a52f7..b6ec7d7c 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m @@ -30,7 +30,7 @@ phip = 1 % sphericity epsilon = 0.451335 % void fraction Ustart = 0.002 Uend = 0.02 -timeStepSize = 0.001; % time interval of pressure data +timeStepSize = 0.0005; % time interval of pressure data Tstart = 0; Tend = t_sim(length(t_sim)); deltaU=(Uend-Ustart)/((Tend-Tstart)/timeStepSize); diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict index 14f3feb9..55a2adba 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict @@ -25,7 +25,7 @@ stopAt endTime; endTime 0.1; -deltaT 0.001; +deltaT 0.0005; writeControl adjustableRunTime; diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/fvSolution b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/fvSolution index 46c6e11d..d6fb7a2b 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/fvSolution +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/fvSolution @@ -41,7 +41,7 @@ solvers relTol 0; } - "(voidfraction|Us|Ksl|dSmoothing)" + "(voidfraction|Us|Ksl|dSmoothing|UsNext|voidfractionNext)" { solver PCG; preconditioner DIC; diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI/DEM/in.liggghts_resume b/tutorials/cfdemSolverPiso/ErgunTestMPI/DEM/in.liggghts_resume index 4ef8baba..47fd1395 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI/DEM/in.liggghts_resume +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI/DEM/in.liggghts_resume @@ -65,7 +65,7 @@ compute_modify thermo_temp dynamic yes #insert the first particles so that dump is not empty dump myDump all stl 1 post/dump_*.stl #run 1 -dump dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius +dump dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius undump myDump run 1 diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties index 2dd8f3d3..558e491e 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties @@ -130,6 +130,7 @@ DiFeliceDragProps velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; + granVelFieldName "Us"; } KochHillDragProps diff --git a/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties index e3d4658f..50e6c9db 100755 --- a/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties @@ -93,6 +93,8 @@ DiFeliceDragProps velFieldName "U"; densityFieldName "rho"; voidfractionFieldName "voidfraction"; + granVelFieldName "Us"; + verbose; } SchillerNaumannDragProps diff --git a/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/system/fvSolution b/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/system/fvSolution index 171310e1..96a1d5c2 100644 --- a/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/system/fvSolution +++ b/tutorials/cfdemSolverPiso/settlingTestMPI/CFD/system/fvSolution @@ -33,7 +33,7 @@ solvers relTol 0; } - U + "(U|k|epsilon|R|nuTilda)" { solver PBiCG; preconditioner DILU; @@ -41,55 +41,13 @@ solvers relTol 0; } - k - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - epsilon - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - R - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - nuTilda - { - solver PBiCG; - preconditioner DILU; - tolerance 1e-05; - relTol 0; - } - - voidfraction + "(voidfraction|Ksl|UsNext|voidfractionNext)" { solver PCG; preconditioner DIC; tolerance 1e-09; relTol 1e-06; } - - Ksl - { - $voidfraction - } - - Us - { - $voidfraction - } } PISO