diff --git a/applications/solvers/cfdemSolverIB/Make/options b/applications/solvers/cfdemSolverIB/Make/options index 9497badb..90f9e40c 100755 --- a/applications/solvers/cfdemSolverIB/Make/options +++ b/applications/solvers/cfdemSolverIB/Make/options @@ -9,12 +9,12 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ -Wno-deprecated-copy EXE_LIBS = \ @@ -27,6 +27,7 @@ EXE_LIBS = \ -ldynamicFvMesh \ -ldynamicMesh \ -lfvOptions \ + -lsampling \ -l$(CFDEM_LIB_NAME) \ $(CFDEM_ADD_LIB_PATHS) \ $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/cfdemSolverIB/cfdemSolverIB.C b/applications/solvers/cfdemSolverIB/cfdemSolverIB.C index cb3f3235..ff5b5d29 100755 --- a/applications/solvers/cfdemSolverIB/cfdemSolverIB.C +++ b/applications/solvers/cfdemSolverIB/cfdemSolverIB.C @@ -6,7 +6,8 @@ Christoph Goniva, christoph.goniva@cfdem.com Copyright (C) 1991-2009 OpenCFD Ltd. Copyright (C) 2009-2012 JKU, Linz - Copyright (C) 2012- DCS Computing GmbH,Linz + Copyright (C) 2012-2015 DCS Computing GmbH,Linz + Copyright (C) 2015- JKU, Linz ------------------------------------------------------------------------------- License This file is part of CFDEMcoupling. @@ -29,11 +30,14 @@ Application Description Transient solver for incompressible flow. - The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, + The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, where additional functionality for CFD-DEM coupling using immersed body (fictitious domain) method is added. Contributions Alice Hager + Daniel Queteschiner + Thomas Lichtenegger + Achuth N. Balachandran Nair \*---------------------------------------------------------------------------*/ @@ -53,23 +57,21 @@ Contributions #include "cellSet.H" +#include "fvOptions.H" // added the fvOptions library + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { + #include "setRootCase.H" - #include "createTime.H" - #include "createDynamicFvMesh.H" - #include "createControl.H" - #include "createTimeControls.H" - #include "createFields.H" - #include "initContinuityErrs.H" + #include "createFvOptions.H" // create cfdemCloud #include "readGravitationalAcceleration.H" @@ -93,24 +95,31 @@ int main(int argc, char *argv[]) // do particle stuff Info << "- evolve()" << endl; - particleCloud.evolve(); + particleCloud.evolve(Us); // Pressure-velocity PISO corrector { + MRF.correctBoundaryVelocity(U); + // Momentum predictor fvVectorMatrix UEqn ( - fvm::ddt(voidfraction,U) + fvm::ddt(voidfraction,U) + MRF.DDt(U) + fvm::div(phi, U) + turbulence->divDevReff(U) + == + fvOptions(U) ); UEqn.relax(); + fvOptions.constrain(UEqn); + if (piso.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); + fvOptions.correct(U); } // --- PISO loop @@ -126,6 +135,7 @@ int main(int argc, char *argv[]) adjustPhi(phi, U, p); + while (piso.correctNonOrthogonal()) { // Pressure corrector @@ -152,12 +162,15 @@ int main(int argc, char *argv[]) } } + laminarTransport.correct(); turbulence->correct(); Info << "particleCloud.calcVelocityCorrection() " << endl; volScalarField voidfractionNext=mesh.lookupObject("voidfractionNext"); particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext); + fvOptions.correct(U); + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/applications/solvers/cfdemSolverIB/createFields.H b/applications/solvers/cfdemSolverIB/createFields.H index 017d4358..ea7e6a25 100755 --- a/applications/solvers/cfdemSolverIB/createFields.H +++ b/applications/solvers/cfdemSolverIB/createFields.H @@ -26,21 +26,6 @@ ), mesh ); - //mod by alice - Info<< "Reading physical velocity field U" << endl; - Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl; - volVectorField Us - ( - IOobject - ( - "Us", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); //======================== // drag law modelling @@ -76,9 +61,8 @@ mesh ); - //mod by alice - Info<< "Reading field phiIB\n" << endl; + Info<< "Reading field voidfraction\n" << endl; volScalarField voidfraction ( IOobject @@ -91,6 +75,21 @@ ), mesh ); + + Info<< "Reading particle velocity field Us\n" << endl; + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + //======================== # include "createPhi.H" @@ -126,3 +125,5 @@ ); //=========================== + +#include "createMRF.H" diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/Make/files b/applications/solvers/cfdemSolverIBContinuousForcing/Make/files new file mode 100755 index 00000000..2635b436 --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/Make/files @@ -0,0 +1,3 @@ +cfdemSolverIBContinuousForcing.C + +EXE=$(CFDEM_APP_DIR)/cfdemSolverIBContinuousForcing diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/Make/options b/applications/solvers/cfdemSolverIBContinuousForcing/Make/options new file mode 100755 index 00000000..90f9e40c --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/Make/options @@ -0,0 +1,34 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +EXE_INC = \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -Wno-deprecated-copy + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools \ + -ldynamicFvMesh \ + -ldynamicMesh \ + -lfvOptions \ + -lsampling \ + -l$(CFDEM_LIB_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) + diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H b/applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H new file mode 100644 index 00000000..a20a5986 --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H @@ -0,0 +1,19 @@ +fvVectorMatrix UEqn +( + fvm::ddt(voidfractionNext,U) + MRF.DDt(U) + + fvm::div(phi, U) + + turbulence->divDevReff(U) + == + fvOptions(U) + + (lambda*(1-voidfractionNext)/U.mesh().time().deltaT())*(fvc::Sp(1,Us)-fvm::Sp(1,U)) +); + +UEqn.relax(); + +fvOptions.constrain(UEqn); + +if (piso.momentumPredictor()) +{ + solve(UEqn == -fvc::grad(p)); + fvOptions.correct(U); +} diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C b/applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C new file mode 100755 index 00000000..4896716f --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C @@ -0,0 +1,129 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Copyright (C) 1991-2009 OpenCFD Ltd. + Copyright (C) 2009-2012 JKU, Linz + Copyright (C) 2012-2015 DCS Computing GmbH,Linz + Copyright (C) 2015- JKU, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with CFDEMcoupling. If not, see . + +Application + cfdemSolverIBContinuousForcing + +Description + Transient solver for incompressible flow. + The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, + where additional functionality for CFD-DEM coupling using immersed body + (fictitious domain) method and a continuous forcing approach is added. +Contributions + Alice Hager + Achuth N. Balachandran Nair +\*---------------------------------------------------------------------------*/ + + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "pisoControl.H" + +#include "cfdemCloudIB.H" +#include "implicitCouple.H" + +#include "averagingModel.H" +#include "regionModel.H" +#include "voidFractionModel.H" + +#include "dynamicFvMesh.H" + +#include "cellSet.H" + +#include "fvOptions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "createFields.H" + #include "initContinuityErrs.H" + #include "createFvOptions.H" + + // create cfdemCloud + #include "readGravitationalAcceleration.H" + cfdemCloudIB particleCloud(mesh); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + //=== dyM =================== + interFace = mag(mesh.lookupObject("voidfractionNext")); + mesh.update(); //dyM + + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" + + // do particle stuff + Info << "- evolve()" << endl; + particleCloud.evolve(Us); + + volScalarField voidfractionNext=mesh.lookupObject("voidfractionNext"); + + // Pressure-velocity PISO corrector + { + MRF.correctBoundaryVelocity(U); + + // Momentum predictor + #include "UEqn.H" + + // --- PISO loop + while (piso.correct()) + { + #include "pEqn.H" + } + } + + laminarTransport.correct(); + turbulence->correct(); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/createFields.H b/applications/solvers/cfdemSolverIBContinuousForcing/createFields.H new file mode 100755 index 00000000..5b15ad97 --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/createFields.H @@ -0,0 +1,143 @@ + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading physical velocity field U" << endl; + Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading particle velocity field Us\n" << endl; + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading the penalization factor field lambda\n" << endl; + volScalarField lambda + ( + IOobject + ( + "lambda", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//======================== +// drag law modelling +//======================== + + Info<< "\nCreating dummy density field rho = 1\n" << endl; + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0) + ); + + + Info<< "Reading field phiIB\n" << endl; + volScalarField phiIB + ( + IOobject + ( + "phiIB", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + //mod by alice + Info<< "Reading field voidfraction\n" << endl; + volScalarField voidfraction + ( + IOobject + ( + "voidfraction", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//======================== + +# include "createPhi.H" + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); + + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + +//=== dyM =================== + + Info<< "Reading field interFace\n" << endl; + volScalarField interFace + ( + IOobject + ( + "interFace", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + //dimensionedScalar("0", dimensionSet(0, -1, 0, 0, 0), 0.0) + dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 0.0) + ); + +//=========================== + +#include "createMRF.H" diff --git a/applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H b/applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H new file mode 100644 index 00000000..4875a1c8 --- /dev/null +++ b/applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H @@ -0,0 +1,35 @@ +volScalarField rUA = 1.0/UEqn.A(); +surfaceScalarField rUAf(fvc::interpolate(rUA)); + +U = rUA*UEqn.H(); + +phi = (fvc::interpolate(U) & mesh.Sf()) + + rUAf*fvc::ddtCorr(U, phi); // Is there additional flux term due to the particle presence? + +MRF.makeRelative(phi); + +adjustPhi(phi, U, p); + +while (piso.correctNonOrthogonal()) +{ + // Pressure corrector + + fvScalarMatrix pEqn + ( + fvm::laplacian(rUA, p) == fvc::div(phi) + particleCloud.ddtVoidfraction() + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(mesh.solver(p.select(piso.finalInnerIter()))); + + if (piso.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } +} + +#include "continuityErrs.H" +U -= rUA*fvc::grad(p); // should we add a pressure correction? +U.correctBoundaryConditions(); +fvOptions.correct(U); diff --git a/doc/CFDEMcoupling_solvers.txt b/doc/CFDEMcoupling_solvers.txt index 49372188..f9ca5779 100644 --- a/doc/CFDEMcoupling_solvers.txt +++ b/doc/CFDEMcoupling_solvers.txt @@ -10,6 +10,7 @@ This section lists all CFDEMcoupling solvers alphabetically. "cfdemSolverIB"_cfdemSolverIB.html, +"cfdemSolverIBContinuousForcing"_cfdemSolverIBContinuousForcing.html, "cfdemSolverMultiphase"_cfdemSolverMultiphase.html, "cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html, "cfdemSolverPiso"_cfdemSolverPiso.html, diff --git a/doc/cfdemSolverIBContinuousForcing.txt b/doc/cfdemSolverIBContinuousForcing.txt new file mode 100644 index 00000000..3db21fbf --- /dev/null +++ b/doc/cfdemSolverIBContinuousForcing.txt @@ -0,0 +1,82 @@ +"CFDEMproject Website"_lws - "Main Page"_main :c + +:link(lws,http://www.cfdem.com) +:link(main,CFDEMcoupling_Manual.html) + +:line + +cfdemSolverIBContinuousForcing command :h3 + +[Description:] + +"cfdemSolverIBContinuousForcing" is a coupled CFD-DEM solver using CFDEMcoupling, +an open-source parallel coupled CFD-DEM framework, for calculating the dynamics +between immersed bodies and the surrounding fluid. Being an implementation of a +continuous forcing immersed boundary method it allows tackling problems where +the body diameter exceeds the maximal size of a fluid cell. + + +Using the toolbox of OpenFOAM®(*) the governing equations of the fluid are +computed and the corrections of velocity and pressure field with respect to the +body-movement information, gained from LIGGGHTS, are incorporated. + + + + +The code of this solver was contributed by A.N. Balachandran Nair, JKU. For more +details, see "Balachandran Nair et al. (2021)"_#BalachandranNair2021 + +[Use:] + +The solver is realized within the open-source framework CFDEMcoupling. Just as +for the unresolved CFD-DEM solver cfdemSolverPiso the file +CFD/constant/couplingProperties contains information about the settings for the +different models. While IOmodel, DataExchangeModel etc. are applicable for all +CFDEMcoupling-solvers, special locate-, force- and void fraction models were +designed for this solver: + +"engineSearchIB"_locateModel_engineSearchIB.html, +"ArchimedesIB"_forceModel_ArchimedesIB.html, +"ShirgaonkarIB"_forceModel_ShirgaonkarIB.html, +"IBVoidfraction"_voidFractionModel_IBVoidFraction.html + +:line + +:link(BalachandranNair2021) +[(Balachandran Nair, 2021)] Balachandran Nair, A.N., Pirker, S. and Saeedipour, M., +"Resolved CFD-DEM simulation of blood flow with a reduced-order RBC model", +Comp. Part. Mech. (2021) + +:line + + +NOTE: +(*) This offering is not approved or endorsed by OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com, and owner of the +OPENFOAM® and OpenCFD® trade marks. +OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and +distributor of the OpenFOAM software via www.openfoam.com. + + + + diff --git a/doc/forceModel_ShirgaonkarIB.txt b/doc/forceModel_ShirgaonkarIB.txt index 313b18c3..63d03b27 100644 --- a/doc/forceModel_ShirgaonkarIB.txt +++ b/doc/forceModel_ShirgaonkarIB.txt @@ -21,12 +21,14 @@ ShirgaonkarIBProps velFieldName "U"; pressureFieldName "pressure"; twoDimensional; + useTorque; treatForceExplicit switch1; \} :pre {U} = name of the finite volume fluid velocity field :ulb,l {pressure} = name of the finite volume pressure field :l {twoDimensional} = optional keyword for conducting a two dimensional calculation :l +{useTorque} = optional keyword for calculating particle torque :l {switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l :ule diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 9d3bb09d..4bffb078 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -10,6 +10,7 @@ cfdemSolverPiso/dir cfdemSolverRhoPimple/dir cfdemSolverIB/dir cfdemSolverPisoScalar/dir +cfdemSolverIBContinuousForcing/dir cfdemSolverRhoPimpleChem/dir cfdemSolverMultiphase/dir cfdemSolverMultiphaseScalar/dir diff --git a/etc/tutorial-list.txt b/etc/tutorial-list.txt index d893faf1..d96302f6 100644 --- a/etc/tutorial-list.txt +++ b/etc/tutorial-list.txt @@ -19,3 +19,5 @@ cfdemSolverIB/twoSpheresGlowinskiMPI/dir cfdemSolverPisoScalar/packedBedTemp/dir cfdemSolverPiso/ErgunTestCG/dir + +cfdemSolverIB/falling_sphere_two_way_coupling/dir diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C index 8fec93d2..09ba105c 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C @@ -36,7 +36,6 @@ Description #include "locateModel.H" #include "dataExchangeModel.H" #include "IOModel.H" -#include #include "IOmanip.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -58,6 +57,7 @@ cfdemCloudIB::cfdemCloudIB angularVelocities_(NULL), pRefCell_(readLabel(mesh.solutionDict().subDict("PISO").lookup("pRefCell"))), pRefValue_(readScalar(mesh.solutionDict().subDict("PISO").lookup("pRefValue"))), + DEMTorques_(NULL), haveEvolvedOnce_(false), skipLagrangeToEulerMapping_(false) { @@ -75,6 +75,7 @@ cfdemCloudIB::cfdemCloudIB cfdemCloudIB::~cfdemCloudIB() { dataExchangeM().destroy(angularVelocities_,3); + dataExchangeM().destroy(DEMTorques_,3); } @@ -82,8 +83,7 @@ cfdemCloudIB::~cfdemCloudIB() void cfdemCloudIB::getDEMdata() { cfdemCloud::getDEMdata(); - Info << "=== cfdemCloudIB::getDEMdata() === particle rotation not considered in CFD" << endl; - //dataExchangeM().getData("omega","vector-atom",angularVelocities_); + dataExchangeM().getData("omega","vector-atom",angularVelocities_); } bool cfdemCloudIB::reAllocArrays() @@ -92,12 +92,24 @@ bool cfdemCloudIB::reAllocArrays() { // get arrays of new length dataExchangeM().allocateArray(angularVelocities_,0.,3); + dataExchangeM().allocateArray(DEMTorques_,0.,3); return true; } return false; } -bool cfdemCloudIB::evolve() +void cfdemCloudIB::giveDEMdata() +{ + cfdemCloud::giveDEMdata(); + dataExchangeM().giveData("hdtorque","vector-atom",DEMTorques_); +} + +inline double ** cfdemCloudIB::DEMTorques() const +{ + return DEMTorques_; +} + +bool cfdemCloudIB::evolve(volVectorField& Us) { numberOfParticlesChanged_ = false; arraysReallocated_=false; @@ -140,6 +152,10 @@ bool cfdemCloudIB::evolve() for (int i=0;i=0) + { + // calc particle velocity + for (int i=0;i<3;i++) rVec[i]=Us.mesh().C()[cell][i]-position(index)[i]; + for (int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i]; + velRot=angVel^rVec; + for (int i=0;i<3;i++) uP[i] = velocities()[index][i]+velRot[i]; + Us[cell] = (1-voidfractions_[index][subCell])*uP; + } + } + } + Us.correctBoundaryConditions(); +} + void cfdemCloudIB::calcVelocityCorrection ( volScalarField& p, @@ -190,7 +234,7 @@ void cfdemCloudIB::calcVelocityCorrection for(int i=0;i<3;i++) uParticle[i] = velocities()[index][i]+velRot[i]; // impose field velocity - U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI]; + U[cellI]= (1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI]; } } //} @@ -213,7 +257,7 @@ void cfdemCloudIB::calcVelocityCorrection U.correctBoundaryConditions(); // correct the pressure as well - p=p+phiIB/U.mesh().time().deltaT(); // do we have to account for rho here? + p=p+phiIB/U.mesh().time().deltaT(); // no need to account for rho since p=(p/rho) in OF p.correctBoundaryConditions(); if (couplingProperties_.found("checkinterface")) diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H index 88041fc8..99439d44 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H @@ -62,6 +62,8 @@ protected: label pRefCell_; scalar pRefValue_; + double **DEMTorques_; + bool haveEvolvedOnce_; bool skipLagrangeToEulerMapping_; @@ -85,7 +87,13 @@ public: bool reAllocArrays(); - bool evolve(); + void giveDEMdata(); + + inline double ** DEMTorques() const; + + bool evolve(volVectorField&); + + void calcForcingTerm(volVectorField&); void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model @@ -96,7 +104,6 @@ public: { return angularVelocities_; } - }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C index a2bbcb56..62b576e3 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C @@ -5,7 +5,8 @@ www.cfdem.com Christoph Goniva, christoph.goniva@cfdem.com Copyright 2009-2012 JKU Linz - Copyright 2012- DCS Computing GmbH, Linz + Copyright 2012-2015 DCS Computing GmbH, Linz + Copyright 2015- JKU Linz ------------------------------------------------------------------------------- License This file is part of CFDEMcoupling. @@ -66,10 +67,14 @@ ShirgaonkarIB::ShirgaonkarIB verbose_(propsDict_.found("verbose")), twoDimensional_(propsDict_.found("twoDimensional")), depth_(1), + multisphere_(propsDict_.found("multisphere")), // drag for a multisphere particle + useTorque_(propsDict_.found("useTorque")), velFieldName_(propsDict_.lookup("velFieldName")), U_(sm.mesh().lookupObject (velFieldName_)), pressureFieldName_(propsDict_.lookup("pressureFieldName")), - p_(sm.mesh().lookupObject (pressureFieldName_)) + p_(sm.mesh().lookupObject (pressureFieldName_)), + solidVolFractionName_(multisphere_?propsDict_.lookup("solidVolFractionName"):word::null), + lambda_(multisphere_?sm.mesh().lookupObject(solidVolFractionName_):volScalarField::null()) { //Append the field names to be probed particleCloud_.probeM().initialize(typeName, typeName+".logDat"); @@ -109,6 +114,7 @@ void ShirgaonkarIB::setForce() const label cellI; vector drag; + vector torque; volVectorField h=forceSubM(0).IBDragPerV(U_,p_); @@ -119,6 +125,7 @@ void ShirgaonkarIB::setForce() const //if(mask[index][0]) //{ drag=vector::zero; + torque=vector::zero; for(int subCell=0;subCell -1) // particle Found { drag += h[cellI]*h.mesh().V()[cellI]; + if(useTorque_) + { + vector rc = particleCloud_.mesh().C()[cellI]; + vector positionCenter = particleCloud_.position(index); + torque += (rc - positionCenter)^h[cellI]*h.mesh().V()[cellI]; + } } } @@ -150,10 +163,38 @@ void ShirgaonkarIB::setForce() const << "," << particleCloud_.impForces()[index][1] << "," << particleCloud_.impForces()[index][2] << endl; + + if(useTorque_) + { + for(int j=0;j<3;j++) particleCloud_.DEMTorques()[index][j] = torque[j]; // Adding to the particle torque; + if(verbose_) Info << "DEMTorques = " << particleCloud_.DEMTorques()[index][0] << "," << particleCloud_.DEMTorques()[index][1] << "," << particleCloud_.DEMTorques()[index][2] << endl; + } //} } + + // Calculate the force if the particle is multisphere template + if(multisphere_) calcForce(); } +void ShirgaonkarIB::calcForce() const +{ + vector dragMS; + + volVectorField h=forceSubM(0).IBDragPerV(U_,p_); + + dragMS = vector::zero; + + forAll(lambda_,cellI) + { + if(lambda_[cellI] > 0) dragMS += h[cellI]*h.mesh().V()[cellI]; + else dragMS = dragMS; + } + + Pout << "Drag force on particle clump = " << dragMS[0] << ", " << dragMS[1] << ", " << dragMS[2] << endl; +} + + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H index a5fb0ec8..b9d4ea8a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.H @@ -5,7 +5,8 @@ www.cfdem.com Christoph Goniva, christoph.goniva@cfdem.com Copyright 2009-2012 JKU Linz - Copyright 2012- DCS Computing GmbH, Linz + Copyright 2012-2015 DCS Computing GmbH, Linz + Copyright 2015- JKU Linz ------------------------------------------------------------------------------- License This file is part of CFDEMcoupling. @@ -66,6 +67,10 @@ private: bool depth_; + const bool multisphere_; + + const bool useTorque_; + word velFieldName_; const volVectorField& U_; @@ -74,6 +79,12 @@ private: const volScalarField& p_; + word solidVolFractionName_; + + const volScalarField& lambda_; + + void calcForce() const; + public: //- Runtime type information diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allclean.sh b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allclean.sh new file mode 100755 index 00000000..3360c516 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allclean.sh @@ -0,0 +1,27 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allclean script for falling sphere test case +# clean CFD and DEM part +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source environment variables +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#- clean up case +echo "deleting data at: $casePath" +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanCase +rm -r $casePath/CFD/clockData +rm $casePath/DEM/post/*.* +#- preserve post directory +touch $casePath/DEM/post/.gitignore +echo "done" + diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allrun.sh b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allrun.sh new file mode 100755 index 00000000..d1b1e2d6 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allrun.sh @@ -0,0 +1,25 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# allrun script for falling sphere testcase +# run falling_sphere_two_way_coupling +# Achuth N. Balachandran Nair - Oct. 2018 +#------------------------------------------------------------------------------ + +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +echo $casePath + +# check if mesh was built +if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then + echo "mesh was built before - using old mesh" +else + echo "mesh needs to be built" + cd CFD/ + blockMesh +fi + +#- run parallel CFD-DEM +bash $casePath/parCFDDEMrun.sh diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/U b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/U new file mode 100755 index 00000000..46fb54a3 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/U @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value $internalField; + } + + channelWall + { + type zeroGradient; + } + + outlet + { + type fixedValue; + value $internalField; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/Us b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/Us new file mode 100755 index 00000000..7f20cd03 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/Us @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object Us; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type zeroGradient; + } + + outlet + { + type zeroGradient; + } + + channelWall + { + type zeroGradient; + } + + frontAndBack + { + type zeroGradient; + } +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/p b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/p new file mode 100755 index 00000000..b4e5f4ed --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/p @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type fixedValue; + value $internalField; + } + + channelWall + { + type zeroGradient; + } + + outlet + { + type fixedValue; + value $internalField; + } + + frontAndBack + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/phiIB b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/phiIB new file mode 100755 index 00000000..307f404f --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/phiIB @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object phiIB; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type fixedValue; + value uniform 0; + } + + channelWall + { + type zeroGradient; + } + + outlet + { + type fixedValue; + value uniform 0; + } + + frontAndBack + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho new file mode 100644 index 00000000..bdd31854 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object rho; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -3 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + inlet + { + type zeroGradient; + } + + channelWall + { + type zeroGradient; + } + + outlet + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/voidfraction b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/voidfraction new file mode 100755 index 00000000..9722d600 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/voidfraction @@ -0,0 +1,44 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object voidfraction; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + inlet + { + type zeroGradient; + } + + outlet + { + type zeroGradient; + } + + channelWall + { + type zeroGradient; + } + + frontAndBack + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/couplingProperties b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/couplingProperties new file mode 100644 index 00000000..94c261e3 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/couplingProperties @@ -0,0 +1,101 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object couplingProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// sub-models & settings + +modelType none; + +checkPeriodicCells; + +verbose; + +couplingInterval 10; + +depth 0; + +voidFractionModel IB; + +locateModel engineIB; + +meshMotionModel noMeshMotion; + +regionModel allRegion; + +dataExchangeModel twoWayMPI; + +IOModel basicIO; + +probeModel off; + +averagingModel dilute; + +clockModel off; + +smoothingModel off; + +forceModels +( + ShirgaonkarIB + ArchimedesIB +); + +momCoupleModels +( +); + +turbulenceModelType "turbulenceProperties"; + +// sub-model properties + +ShirgaonkarIBProps +{ + velFieldName "U"; + pressureFieldName "p"; + densityFieldName "rho"; + verbose; + useTorque; +} + +ArchimedesIBProps +{ + densityFieldName "rho"; + gravityFieldName "g"; + voidfractionFieldName "voidfractionNext"; +} + +twoWayMPIProps +{ + liggghtsPath "../DEM/in.liggghts_run"; +} + +IBProps +{ + maxCellsPerParticle 20000; + alphaMin 0.30; + scaleUpVol 1.0; +} + +engineIBProps +{ + engineProps + { + treeSearch false; + } + zSplit 8; + xySplit 16; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/dynamicMeshDict b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/dynamicMeshDict new file mode 100644 index 00000000..017c6dfd --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/dynamicMeshDict @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object dynamicMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dynamicFvMesh staticFvMesh; + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/g b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/g new file mode 100755 index 00000000..09cf6c45 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 -9.81 0); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/liggghtsCommands new file mode 100755 index 00000000..cf03eae8 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/liggghtsCommands @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object liggghtsCommands; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +liggghtsCommandModels +( + runLiggghts +); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/transportProperties b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/transportProperties new file mode 100755 index 00000000..d22bf944 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/transportProperties @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +transportModel Newtonian; + +// Fluid properties of plasma // dextran +nu nu [ 0 2 -1 0 0 0 0 ] 5e-06; + +CrossPowerLawCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + m m [ 0 0 1 0 0 0 0 ] 1; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +BirdCarreauCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [ 0 0 1 0 0 0 0 ] 0; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/turbulenceProperties b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/turbulenceProperties new file mode 100755 index 00000000..6d46e885 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/turbulenceProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/octave/plot_data.m b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/octave/plot_data.m new file mode 100644 index 00000000..660a7b61 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/octave/plot_data.m @@ -0,0 +1,33 @@ +clear all; +clc; + +%% Script to plot the angular velocities and the particle positions + +A = importdata('../../DEM/post/angular_velocity_no_coupling.txt',' ',1); +B = importdata('../../DEM/post/position_no_coupling.txt',' ',1); +C = importdata('../../DEM/post/angular_velocity.txt',' ',1); +D = importdata('../../DEM/post/position.txt',' ',1); + +pos1 = B.data(); +omega1 = A.data(); +pos2 = D.data(); +omega2 = C.data(); + +time = omega1(:,1); +omegax1 = omega1(:,2); +omegay1 = omega1(:,3); +omegaz1 = omega1(:,4); + +time2 = omega2(:,1); +omegax2 = omega2(:,2); +omegay2 = omega2(:,3); +omegaz2 = omega2(:,4); + +figure +plot(time,omegaz1,'-.-',time2,omegaz2,'Linewidth',1.5) +xlabel('Time (s)') +ylabel('Angular velocity (1/s)') +legend('One-way coupling','Two-way coupling') +axis([0 0.5 0 11]) +set(gca,'FontSize',12) +print('angular_velocity_compare.eps') diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/blockMeshDict b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/blockMeshDict new file mode 100644 index 00000000..12997756 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/blockMeshDict @@ -0,0 +1,77 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + (0 0 0) + (0.025 0 0) + (0.025 0.1 0) + (0 0.1 0) + (0 0 0.025) + (0.025 0 0.025) + (0.025 0.1 0.025) + (0 0.1 0.025) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (20 80 20) simpleGrading (1 1 1) +); + +edges +( +); + +boundary +( + inlet + { + type wall; + faces + ( + (0 3 2 1) + ); + } + + channelWall + { + type wall; + faces + ( + (2 6 5 1) + (0 4 7 3) + (3 7 6 2) + (1 5 4 0) + ); + } + + outlet + { + type wall; + faces + ( + (4 5 6 7) + ); + } +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/controlDict b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/controlDict new file mode 100644 index 00000000..a8cfa722 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/controlDict @@ -0,0 +1,56 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application cfdemSolverIB; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 0.5; + +deltaT 0.0005; + +writeControl adjustableRunTime; + +writeInterval 0.005; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable off; + +adjustTimeStep no; + +maxCo 0.6; + +maxDeltaT 1; + + +// ************************************************************************* // + diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/decomposeParDict b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/decomposeParDict new file mode 100644 index 00000000..aa2a0618 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/decomposeParDict @@ -0,0 +1,32 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + note "mesh decomposition control dictionary"; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 4; + +//- Keep owner and neighbour on same processor for faces in zones: +// preserveFaceZones (heater solid1 solid3); + + method simple; + +simpleCoeffs +{ + n (1 1 4); + delta 0.001; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvOptions b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvOptions new file mode 100644 index 00000000..1b83175a --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvOptions @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object fvOptions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +momentumSource +{ + type meanVelocityForce; + active yes; + + meanVelocityForceCoeffs + { + selectionMode all; + + fields (U); + Ubar (0.01 0 0); + relaxation 0.2; + } +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSchemes b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSchemes new file mode 100755 index 00000000..433500f0 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSchemes @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; + grad(p) Gauss linear; + grad(U) Gauss linear; +} + +divSchemes +{ + default Gauss linear; + div(phi,U) Gauss limitedLinearV 1; + div(phi,k) Gauss limitedLinear 1; + div(phi,epsilon) Gauss limitedLinear 1; + div(phi,R) Gauss limitedLinear 1; + div(R) Gauss linear; + div(phi,nuTilda) Gauss limitedLinear 1; + div((nuEff*dev(grad(U).T()))) Gauss linear; + div(U) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(nuEff,U) Gauss linear corrected; + laplacian((1|A(U)),p) Gauss linear corrected; + laplacian((voidfraction2|A(U)),p) Gauss linear corrected; + laplacian(DkEff,k) Gauss linear corrected; + laplacian(DepsilonEff,epsilon) Gauss linear corrected; + laplacian(DREff,R) Gauss linear corrected; + laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; + laplacian(phiIB) Gauss linear corrected; + laplacian(U) Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; + interpolate(U) linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSolution b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSolution new file mode 100755 index 00000000..ebf44623 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSolution @@ -0,0 +1,110 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + pFinal + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + U + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + 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; + } + + phiIB + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + "(U|k)Final" + { + $U; + tolerance 1e-05; + relTol 0; + } + +} + +PISO +{ + nCorrectors 4; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + +PIMPLE +{ + nOuterCorrectors 3; + nCorrectors 2; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/in.liggghts_run b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/in.liggghts_run new file mode 100644 index 00000000..004d3794 --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/in.liggghts_run @@ -0,0 +1,88 @@ + +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +variable deltat equal 5e-05 + +atom_style sphere +atom_modify map array sort 0 0 +boundary m m m +newton off + +communicate single vel yes +processors 1 4 1 + +units si + +region reg block 0 0.025 0 0.1 0 0.025 units box +create_box 1 reg + +neighbor 0.1 bin +neigh_modify delay 0 binsize 0.01 + +# material properties required for granular pair styles +fix m1 all property/global youngsModulus peratomtype 5e05 +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 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +timestep ${deltat} + +fix gravi all gravity 9.81 vector 0.0 -1.0 0.0 + +# walls +fix xwall1 all wall/gran model hertz tangential history primitive type 1 xplane 0.0 +fix xwall2 all wall/gran model hertz tangential history primitive type 1 xplane 0.025 +fix zwall1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0 +fix zwall2 all wall/gran model hertz tangential history primitive type 1 zplane 0.025 +fix ywall1 all wall/gran model hertz tangential history primitive type 1 yplane 0.0 +fix ywall2 all wall/gran model hertz tangential history primitive type 1 yplane 0.1 + +# cfd coupling +fix cfd all couple/cfd couple_every 10 mpi +fix cfd2 all couple/cfd/force + +# create atom +create_atoms 1 single 0.0125 0.075 0.0125 units box +set atom 1 diameter 0.01 density 1.1 vx 0 vy 0 vz 0 omegaz 10 + + +# getting the angular velocities to print out +variable omega1 equal omegax[1] +variable omega2 equal omegay[1] +variable omega3 equal omegaz[1] + +# getting the particle position to print out +variable x1 equal x[1] +variable y1 equal y[1] +variable z1 equal z[1] + +variable time equal step*dt + +fix extra1 all print 10 "${time} ${x1} ${y1} ${z1}" file ../DEM/post/position.txt screen no +fix extra2 all print 10 "${time} ${omega1} ${omega2} ${omega3}" file ../DEM/post/angular_velocity.txt screen no + +# integrator +fix integr all nve/sphere + +# screen output +compute rke all erotate/sphere +thermo_style custom step atoms ke c_rke vol +thermo 1000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic no extra 0 + +# insert the first particles so that dump is not empty +dump dmp1 all custom 100 ../DEM/post/dump.liggghts_run id type x y z vx vy vz & + omegax omegay omegaz radius fx fy fz + +#force : f_couple_cfd[0] f_couple_cfd[1] f_couple_cfd[2] +#node : f_couple_cfd[6] +#cell id : f_couple_cfd[7] + +run 1 + diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/post/.gitignore b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/post/.gitignore new file mode 100644 index 00000000..e69de29b diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/parCFDDEMrun.sh b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/parCFDDEMrun.sh new file mode 100644 index 00000000..7df8a47b --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/parCFDDEMrun.sh @@ -0,0 +1,74 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# parCFDDEMrun script for falling sphere testcase +# run falling_sphere_two_way_coupling +# Achuth N. Balachandran Nair - Oct. 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- include functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath=$casePath +headerText="run_parallel_CFDEMIB_two_way_coupling" +logfileName="log_$headerText" +solverName="cfdemSolverIB" +nrProcs="4" +machineFileName="none" # yourMachinefileName | none +debugMode="off" # on | off| strict +testHarnessPath="$CFDEM_TEST_HARNESS_PATH" +runOctave="false" +postproc="false" +#------------------------------------------------------------------------------ + +#- call function to run a parallel CFD-DEM case +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode + +if [ $runOctave == "true" ] + then + + cd $casePath/CFD/octave + octave plot_data.m + evince angular_velocity_compare.eps +fi + +if [ $postproc == "true" ] + then + #- get VTK data from liggghts dump file + cd $casePath/DEM/post + lpp dump* + + #- get VTK data from CFD sim + cd $casePath/CFD + reconstructParMesh -constant -mergeTol 1e-06 + reconstructPar -noLagrangian + foamToVTK + + #- start paraview + paraview + + #- keep terminal open (if started in new terminal) + echo "...press enter to clean up case" + echo "press Ctr+C to keep data" + read +fi + +#- copy log file to test harness +cp $casePath/$logfileName $testHarnessPath + +#- clean up case +echo "deleting data at: $casePath" +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanCase +rm -r $casePath/CFD/clockData +rm $casePath/DEM/post/*.* +#- preserve post directory +touch $casePath/DEM/post/.gitignore +echo "done" + diff --git a/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/post.sh b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/post.sh new file mode 100755 index 00000000..dc6b9fce --- /dev/null +++ b/tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/post.sh @@ -0,0 +1,32 @@ +#!/bin/bash + +#===================================================================# +# post run script for testcase +# postprocess falling_sphere_two_way_coupling +# Achuth N. Balachandran Nair - Oct. 2018 +#===================================================================# + +#- source CFDEM env vars +. ~/.bashrc + +#- include functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#--------------------------------------------------------------------------------# +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +#--------------------------------------------------------------------------------# + +#- get VTK data from liggghts dump file +cd $casePath/DEM/post +lpp dump* + +#- get VTK data from CFD sim +cd $casePath/CFD +reconstructParMesh -constant -mergeTol 1e-06 +reconstructPar -noLagrangian +foamToVTK + +#- start paraview +paraview + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allclean.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allclean.sh new file mode 100755 index 00000000..566e5b4a --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allclean.sh @@ -0,0 +1,31 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allclean script for redBloodCellPoiseuilleFlow test case +# clean up redBloodCellPoiseuilleFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#------------------------------------------------------------------------------ +#- clean up case +echo "deleting data at: $casePath" +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanTimeDirectories +rm -rf processor* > /dev/null 2>&1 +rm -rf clockData > /dev/null 2>&1 +echo "Cleaning $casePath/DEM/post" +rm $casePath/DEM/post/*.* > /dev/null 2>&1 +#rm -r $casePath/DEM/post/restart/*.* +#touch $casePath/DEM/post/restart/.gitignore +echo "done" + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allrun.sh new file mode 100755 index 00000000..0c6383b0 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allrun.sh @@ -0,0 +1,30 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allrun script for redBloodCellPoiseuilleFlow test case +# run redBloodCellPoiseuilleFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +# check if mesh was built +if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then + echo "mesh was built before - using old mesh" +else + echo "mesh needs to be built" + cd $casePath/CFD + blockMesh +fi + +if [ -f "$casePath/DEM/post/restart/liggghts.restart" ]; then + echo "LIGGGHTS init was run before - using existing restart file" +else + #- run serial DEM + $casePath/DEMrun.sh +fi + +#- run parallel CFD-DEM +bash $casePath/parCFDDEMrun.sh diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/U b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/U new file mode 100644 index 00000000..fa7ec155 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/U @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); // velocity in micrometres per second + +boundaryField +{ + inlet + { + type cyclicAMI; + } + + outlet + { + type cyclicAMI; + } + + channelWall + { + type noSlip; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/Us b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/Us new file mode 100644 index 00000000..61469e35 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/Us @@ -0,0 +1,42 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object Us; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type cyclicAMI; + value uniform (0 0 0); + } + + outlet + { + type cyclicAMI; + value uniform (0 0 0); + } + + channelWall + { + type noSlip; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/lambda b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/lambda new file mode 100644 index 00000000..145eb79c --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/lambda @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object lambda; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 10; // Value to be determined based on CFD time-step. +// The ratio between lambda/dt > 1 at all times to achieve proper forcing of the momentum equation. + +boundaryField +{ + inlet + { + type cyclicAMI; + } + + outlet + { + type cyclicAMI; + } + + channelWall + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/p b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/p new file mode 100644 index 00000000..258e96c2 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/p @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type cyclicAMI; + } + + outlet + { + type cyclicAMI; + } + + channelWall + { + type zeroGradient; + } + +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/phiIB b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/phiIB new file mode 100644 index 00000000..fdec6722 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/phiIB @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object phiIB; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type cyclicAMI; + value uniform 0; + } + + outlet + { + type cyclicAMI; + value uniform 0; + } + + channelWall + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/voidfraction b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/voidfraction new file mode 100644 index 00000000..31e7f7b1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/voidfraction @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object voidfraction; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + inlet + { + type cyclicAMI; + } + + outlet + { + type cyclicAMI; + } + + channelWall + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/couplingProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/couplingProperties new file mode 100644 index 00000000..aa8dfa5e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/couplingProperties @@ -0,0 +1,91 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object couplingProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// sub-models & settings + +modelType none; + +checkPeriodicCells; + +couplingInterval 10; + +depth 0; + +voidFractionModel IB; + +locateModel engineIB; + +meshMotionModel noMeshMotion; + +regionModel allRegion; + +dataExchangeModel twoWayMPI; + +IOModel off;//basicIO; + +probeModel off; + +averagingModel dilute; + +clockModel off; + +smoothingModel off; + +forceModels +( + ShirgaonkarIB +); + +momCoupleModels +( +); + +turbulenceModelType "turbulenceProperties"; + + +// sub-model properties + +ShirgaonkarIBProps +{ + velFieldName "U"; + pressureFieldName "p"; + densityFieldName "rho"; + solidVolFractionName "lambda"; +} + +IBProps +{ + maxCellsPerParticle 20000; + alphaMin 0.30; + scaleUpVol 1.0; +} + +engineIBProps +{ + engineProps + { + treeSearch true; + } + zSplit 8; + xySplit 16; +} + +twoWayMPIProps +{ + liggghtsPath "../DEM/in.liggghts_run"; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/g b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/g new file mode 100644 index 00000000..d4b9882c --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 0 0); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/liggghtsCommands new file mode 100644 index 00000000..93a4de77 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/liggghtsCommands @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object liggghtsCommands; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +liggghtsCommandModels +( + runLiggghts +); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/transportProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/transportProperties new file mode 100644 index 00000000..c1b4ee34 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/transportProperties @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +transportModel Newtonian; + +// Fluid properties of plasma // dextran +nu nu [ 0 2 -1 0 0 0 0 ] 1.2; //4.9524; + +CrossPowerLawCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + m m [ 0 0 1 0 0 0 0 ] 1; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +BirdCarreauCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [ 0 0 1 0 0 0 0 ] 0; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/turbulenceProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/turbulenceProperties new file mode 100644 index 00000000..f6d9d4cb --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/turbulenceProperties @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/blockMeshDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/blockMeshDict new file mode 100644 index 00000000..d5bf05e7 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/blockMeshDict @@ -0,0 +1,123 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + (-2.65 -2.65 0) + ( 2.65 -2.65 0) + ( 2.65 2.65 0) + (-2.65 2.65 0) + (-3.5355 -3.5355 0) + ( 3.5355 -3.5355 0) + ( 3.5355 3.5355 0) + (-3.5355 3.5355 0) + + (-2.65 -2.65 90) + ( 2.65 -2.65 90) + ( 2.65 2.65 90) + (-2.65 2.65 90) + (-3.5355 -3.5355 90) + ( 3.5355 -3.5355 90) + ( 3.5355 3.5355 90) + (-3.5355 3.5355 90) +); + +blocks +( + // inner block + hex ( 0 1 2 3 8 9 10 11) (8 8 93) simpleGrading (1 1 1) + // slice 1 + hex ( 0 4 5 1 8 12 13 9) (3 8 93) simpleGrading (1 1 1) + // slice 2 + hex ( 1 5 6 2 9 13 14 10) (3 8 93) simpleGrading (1 1 1) + // slice 3 + hex ( 2 6 7 3 10 14 15 11) (3 8 93) simpleGrading (1 1 1) + // slice 4 + hex ( 3 7 4 0 11 15 12 8) (3 8 93) simpleGrading (1 1 1) +); + + + +edges +( + arc 4 5 ( 0 -5 0) + arc 12 13 ( 0 -5 90) + + arc 5 6 ( 5 0 0) + arc 13 14 ( 5 0 90) + + arc 6 7 ( 0 5 0) + arc 14 15 ( 0 5 90) + + arc 7 4 (-5 0 0) + arc 15 12 (-5 0 90) +); + +boundary +( + inlet + { + type cyclicAMI; + neighbourPatch outlet; + faces + ( + (0 3 2 1) + (0 1 5 4) + (1 2 6 5) + (2 3 7 6) + (3 0 4 7) + ); + transform translational; + separationVector ( 0 0 90); + } + + outlet + { + type cyclicAMI; + neighbourPatch inlet; + faces + ( + (8 9 10 11) + (13 9 8 12) + (14 10 9 13) + (10 14 15 11) + (12 8 11 15) + ); + transform translational; + separationVector ( 0 0 -90); + } + + channelWall + { + type wall; + faces + ( + (13 12 4 5) + (5 6 14 13) + (6 7 15 14) + (7 4 12 15) + ); + } +); + +mergePatchPairs +( +); + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/controlDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/controlDict new file mode 100644 index 00000000..ade3b7e2 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/controlDict @@ -0,0 +1,121 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application pisoFoam; + +startFrom latestTime; + +startTime 0; + +stopAt endTime; + +endTime 500000; + +deltaT 10; + +writeControl adjustableRunTime; + +writeInterval 10000; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep no; + +maxCo 1.0; + +maxDeltaT 10; + + +// ************************************************************************* // +DimensionedConstants +{ + unitSet micro; // SI; + + SICoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-11; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-34; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-31; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-27; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-27; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } + + microCoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-20; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-13; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-16; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-12; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-12; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } +} diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/decomposeParDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/decomposeParDict new file mode 100644 index 00000000..1b221387 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/decomposeParDict @@ -0,0 +1,28 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 2; + +method simple; + +simpleCoeffs +{ + n (1 1 2); + delta 0.001; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvOptions b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvOptions new file mode 100644 index 00000000..b743457c --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvOptions @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object fvOptions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +momentumSource +{ + type meanVelocityForce; + active yes; + + meanVelocityForceCoeffs + { + selectionMode all; + fields (U); + // Note: set mean velocity here + Ubar ( 0 0 0.0008 ); + relaxation 0.9; + } +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSchemes b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSchemes new file mode 100644 index 00000000..005e8180 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSchemes @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss skewCorrected linear; + grad(p) Gauss skewCorrected linear; + grad(U) Gauss skewCorrected linear; +} + +divSchemes +{ + default Gauss skewCorrected linear; + div(phi,U) Gauss skewCorrected upwind; + div(phi,k) Gauss limitedLinear 1; + div(phi,epsilon) Gauss limitedLinear 1; + div(phi,R) Gauss limitedLinear 1; + div(R) Gauss linear; + div(phi,nuTilda) Gauss limitedLinear 1; + div((nuEff*dev(grad(U).T()))) Gauss linear; + div(U) Gauss skewCorrected lineaUpwind; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(nuEff,U) Gauss linear corrected; + laplacian((1|A(U)),p) Gauss linear corrected; + laplacian((voidfraction2|A(U)),p) Gauss linear corrected; + laplacian(DkEff,k) Gauss linear corrected; + laplacian(DepsilonEff,epsilon) Gauss linear corrected; + laplacian(DREff,R) Gauss linear corrected; + laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; + laplacian(phiIB) Gauss linear corrected; + laplacian(U) Gauss linear corrected; +} + +interpolationSchemes +{ + default skewCorrected linear; + interpolate(U) skewCorrected linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSolution b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSolution new file mode 100644 index 00000000..7d964cc3 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSolution @@ -0,0 +1,116 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + pFinal + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + U + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + 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; + } + + phiIB + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + "(U|k)Final" + { + $U; + tolerance 1e-05; + relTol 0; + } + + Us + { + $U; + tolerance 1e-05; + relTol 0; + } + +} + +PISO +{ + nCorrectors 2; + nNonOrthogonalCorrectors 1; + pRefCell 0; + pRefValue 0; +} + +relaxationFactors +{ + equations + { + ".*" 0.8; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/data/spheres.txt b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/data/spheres.txt new file mode 100644 index 00000000..55c1fc96 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/data/spheres.txt @@ -0,0 +1,10 @@ +1.58162822 1.14912017 0.00000000 1.95500000 +0.60412822 1.85931549 0.00000000 1.95500000 +-0.60412822 1.85931549 0.00000000 1.95500000 +-1.58162822 1.14912017 0.00000000 1.95500000 +-1.95500000 0.00000000 0.00000000 1.95500000 +-1.58162822 -1.14912017 0.00000000 1.95500000 +-0.60412822 -1.85931549 0.00000000 1.95500000 +0.60412822 -1.85931549 0.00000000 1.95500000 +1.58162822 -1.14912017 0.00000000 1.95500000 +1.95500000 -0.00000000 0.00000000 1.95500000 diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_init b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_init new file mode 100644 index 00000000..ee7ed4f4 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_init @@ -0,0 +1,116 @@ +# Insert a RBC into a cylinder, then induce flow +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 300 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Sben equal ${Sn}*10 # analogous bending modulus of the parallel bond +variable Stor equal ${Sn}*10 # analogous torsional modulus of the parallel bond + + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary f f p + +# communication +newton off +communicate single vel yes + +# units +units micro + +# region and box creation +region domain block -10 10 -10 10 0 90 units box +create_box 2 domain + +# neighbors +neighbor 1.0 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e03 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.5 +fix m3 all property/global coefficientRestitution peratomtypepair 2 1e-06 1e-03 1e-03 1e-06 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.1 0.1 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 2.65 0.0 1 1e50 1e50 + +# place an infinite cylindrical wall with radius 5 µm in the middle of the domain (0,0) +# the cylinder's axis is aligend along the z-axis ('zcylinder') +fix cylwalls all wall/gran model hertz tangential history primitive type 2 zcylinder 5.0 0.0 0.0 + +# particle template, insertion region and bond formation +fix pts1 all particletemplate/multiplespheres 15485863 atom_type 1 & +density constant 1098 nspheres 10 ntry 1000000 spheres file ../DEM/data/spheres.txt scale 1 & +bonded yes/explicit & +nbond_pairs 10 & +1 2 & +2 3 & +3 4 & +4 5 & +5 6 & +6 7 & +7 8 & +8 9 & +9 10 & +10 1 & +bond_type 1 + +fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0 + +region bc1 block -10 10 -10 10 10 70 units box + +fix ins1 all insert/pack seed 32452867 distributiontemplate pdd1 & + maxattempt 10000 insert_every once overlapcheck yes vel constant 0.0 0.0 0.0 & + region bc1 particles_in_region 1 pos 0 0 60 ntry_mc 100000 + +run 1000 + +# set timestep +timestep 0.001 + +# check timestep +fix ts_check all check/timestep/gran 1000 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 10000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump*.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bond type f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds*.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 100000 + +write_restart ../DEM/post/restart/liggghts.restart + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_run b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_run new file mode 100644 index 00000000..223c5fb2 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_run @@ -0,0 +1,133 @@ +# RBC in a cylinder moved by fluid +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 25 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Ab equal (22/7)*${rb}${rb} # bond area +variable Sben equal ${Sn}/${Ab} # analogous bending modulus of the parallel bond +variable Stor equal ${St}/${Ab} # analogous torsional modulus of the parallel bond + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary f f p + +# communication +newton off +communicate single vel yes + +# units and processors +units micro +processors 1 1 2 + +# read the restart file +read_restart ../DEM/post/restart/liggghts.restart + +# neighbors +neighbor 1.0 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e04 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.18 +fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.9 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.3 0.3 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 2.65 2.0 1.0 1e50 1e50 + +# place an infinite cylindrical wall with radius 5 µm in the middle of the domain (0,0) +# the cylinder's axis is aligend along the z-axis ('zcylinder') +fix cylwalls all wall/gran model hertz tangential history primitive type 2 zcylinder 5.0 0.0 0.0 + +# reset timestep after reading the restart file +reset_timestep 0 +timestep 1 + +# cfd coupling +fix cfd all couple/cfd couple_every 10 mpi +fix cfd2 all couple/cfd/force + +# apply nve integration to all constituent spheres - single sphere treatment +fix integr all nve/sphere + +compute unwra all property/atom xu yu zu + +# variables for data processing +variable x1 equal c_unwra[1][1] +variable y1 equal c_unwra[1][2] +variable z1 equal c_unwra[1][3] +variable x2 equal c_unwra[2][1] +variable y2 equal c_unwra[2][2] +variable z2 equal c_unwra[2][3] +variable x3 equal c_unwra[3][1] +variable y3 equal c_unwra[3][2] +variable z3 equal c_unwra[3][3] +variable x4 equal c_unwra[4][1] +variable y4 equal c_unwra[4][2] +variable z4 equal c_unwra[4][3] +variable x5 equal c_unwra[5][1] +variable y5 equal c_unwra[5][2] +variable z5 equal c_unwra[5][3] +variable x6 equal c_unwra[6][1] +variable y6 equal c_unwra[6][2] +variable z6 equal c_unwra[6][3] +variable x7 equal c_unwra[7][1] +variable y7 equal c_unwra[7][2] +variable z7 equal c_unwra[7][3] +variable x8 equal c_unwra[8][1] +variable y8 equal c_unwra[8][2] +variable z8 equal c_unwra[8][3] +variable x9 equal c_unwra[9][1] +variable y9 equal c_unwra[9][2] +variable z9 equal c_unwra[9][3] +variable x10 equal c_unwra[10][1] +variable y10 equal c_unwra[10][2] +variable z10 equal c_unwra[10][3] +variable time equal step*dt + +fix extra1 all print 10 "${time} ${x1} ${y1} ${z1} ${x2} ${y2} ${z2} ${x3} ${y3} ${z3} ${x4} ${y4} ${z4} ${x5} ${y5} ${z5} ${x6} ${y6} ${z6} ${x7} ${y7} ${z7} ${x8} ${y8} ${z8} ${x9} ${y9} ${z9} ${x10} ${y10} ${z10}" & + file ../DEM/post/particle_position.txt screen no + +# check timestep +fix ts_check all check/timestep/gran 1 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 1000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bond type f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 1 + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/post/restart/.gitignore b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/post/restart/.gitignore new file mode 100644 index 00000000..e69de29b diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEMrun.sh new file mode 100755 index 00000000..cdae4cfb --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEMrun.sh @@ -0,0 +1,26 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# DEMrun script for redBloodCellPoiseuilleFlow test case +# init redBloodCellPoiseuilleFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- include functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath="$casePath" +headerText="run_liggghts_init_DEM" +logfileName="log_$headerText" +solverName="in.liggghts_init" +debugMode="off" +#------------------------------------------------------------------------------ + +#- call function to run DEM case +DEMrun $logpath $logfileName $casePath $headerText $solverName $debugMode + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/parCFDDEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/parCFDDEMrun.sh new file mode 100755 index 00000000..2ef4169d --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/parCFDDEMrun.sh @@ -0,0 +1,59 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# parCFDDEMrun script for redBloodCellPoiseuilleFlow test case +# run redBloodCellPoiseuilleFlow CFD-DEM +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath=$casePath +headerText="run_parallel_cfdemSolverIBContinuousForcing_redBloodCellPoiseuilleFlow" +logfileName="log_$headerText" +solverName="cfdemSolverIBContinuousForcing" +nrProcs="2" +machineFileName="none" # yourMachinefileName | none +debugMode="off" # on | off| strict +testHarnessPath="$CFDEM_TEST_HARNESS_PATH" +runPython="false" +postproc="false" +#------------------------------------------------------------------------------ + +#- call function to run a parallel CFD-DEM case +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode + +if [ $runPython == "true" ] + then + cd $casePath + python results.py +fi + +if [ $postproc == "true" ] + then + #- get VTK data from liggghts dump file + cd $casePath/DEM/post + lpp dump* + + #- get VTK data from CFD sim + cd $casePath/CFD + reconstructParMesh -constant -mergeTol 1e-06 + reconstructPar -noLagrangian + foamToVTK + + #- start paraview + paraview + + #- keep terminal open (if started in new terminal) + echo "...press enter to clean up case" + echo "press Ctr+C to keep data" + read +fi + + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/results.py b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/results.py new file mode 100755 index 00000000..c94da0c9 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/results.py @@ -0,0 +1,57 @@ +#!/usr/bin/python + +import numpy as np +import matplotlib.pyplot as plt + +plt.rcParams["font.family"] = "serif" + +# Reading sphere data file +t,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,x9,y9,z9,x10,y10,z10 = np.genfromtxt('DEM/post/particle_position.txt', skip_header=0,usecols = (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30),unpack=True) +t = t*1e-06 + +# Calculating the centre of mass of RBC +M = 154324.770275 +m = (4/3)*np.pi*1.955**3 +xcm = (x1+x2+x3+x4+x5+x6+x7+x8+x9+x10)/10 +ycm = (y1+y2+y3+y4+y5+y6+y7+y8+y9+y10)/10 +zcm = (z1+z2+z3+z4+z5+z6+z7+z8+z9+z10)/10 + +# Calculating the gyration tensor +G11 = ((x1-xcm)*(x1-xcm) + (x2-xcm)*(x2-xcm) + (x3-xcm)*(x3-xcm) + (x4-xcm)*(x4-xcm) + (x5-xcm)*(x5-xcm) + (x6-xcm)*(x6-xcm) + (x7-xcm)*(x7-xcm) + (x8-xcm)*(x8-xcm) + (x9-xcm)*(x9-xcm) + (x10-xcm)*(x10-xcm))/10 +G12 = ((x1-xcm)*(y1-ycm) + (x2-xcm)*(y2-ycm) + (x3-xcm)*(y3-ycm) + (x4-xcm)*(y4-ycm) + (x5-xcm)*(y5-ycm) + (x6-xcm)*(y6-ycm) + (x7-xcm)*(y7-ycm) + (x8-xcm)*(y8-ycm) + (x9-xcm)*(y9-ycm) + (x10-xcm)*(y10-ycm))/10 +G13 = ((x1-xcm)*(z1-zcm) + (x2-xcm)*(z2-zcm) + (x3-xcm)*(z3-zcm) + (x4-xcm)*(z4-zcm) + (x5-xcm)*(z5-zcm) + (x6-xcm)*(z6-zcm) + (x7-xcm)*(z7-zcm) + (x8-xcm)*(z8-zcm) + (x9-xcm)*(z9-zcm) + (x10-xcm)*(z10-zcm))/10 +G21 = ((y1-ycm)*(x1-xcm) + (y2-ycm)*(x2-xcm) + (y3-ycm)*(x3-xcm) + (y4-ycm)*(x4-xcm) + (y5-ycm)*(x5-xcm) + (y6-ycm)*(x6-xcm) + (y7-ycm)*(x7-xcm) + (y8-ycm)*(x8-xcm) + (y9-ycm)*(x9-xcm) + (y10-ycm)*(x10-xcm))/10 +G22 = ((y1-ycm)*(y1-ycm) + (y2-ycm)*(y2-ycm) + (y3-ycm)*(y3-ycm) + (y4-ycm)*(y4-ycm) + (y5-ycm)*(y5-ycm) + (y6-ycm)*(y6-ycm) + (y7-ycm)*(y7-ycm) + (y8-ycm)*(y8-ycm) + (y9-ycm)*(y9-ycm) + (y10-ycm)*(y10-ycm))/10 +G23 = ((y1-ycm)*(z1-zcm) + (y2-ycm)*(z2-zcm) + (y3-ycm)*(z3-zcm) + (y4-ycm)*(z4-zcm) + (y5-ycm)*(z5-zcm) + (y6-ycm)*(z6-zcm) + (y7-ycm)*(z7-zcm) + (y8-ycm)*(z8-zcm) + (y9-ycm)*(z9-zcm) + (y10-ycm)*(z10-zcm))/10 +G31 = ((z1-zcm)*(x1-xcm) + (z2-zcm)*(x2-xcm) + (z3-zcm)*(x3-xcm) + (z4-zcm)*(x4-xcm) + (z5-zcm)*(x5-xcm) + (z6-zcm)*(x6-xcm) + (z7-zcm)*(x7-xcm) + (z8-zcm)*(x8-xcm) + (z9-zcm)*(x9-xcm) + (z10-zcm)*(x10-xcm))/10 +G32 = ((z1-zcm)*(y1-ycm) + (z2-zcm)*(y2-ycm) + (z3-zcm)*(y3-ycm) + (z4-zcm)*(y4-ycm) + (z5-zcm)*(y5-ycm) + (z6-zcm)*(y6-ycm) + (z7-zcm)*(y7-ycm) + (z8-zcm)*(y8-ycm) + (z9-zcm)*(y9-ycm) + (z10-zcm)*(y10-ycm))/10 +G33 = ((z1-zcm)*(z1-zcm) + (z2-zcm)*(z2-zcm) + (z3-zcm)*(z3-zcm) + (z4-zcm)*(z4-zcm) + (z5-zcm)*(z5-zcm) + (z6-zcm)*(z6-zcm) + (z7-zcm)*(z7-zcm) + (z8-zcm)*(z8-zcm) + (z9-zcm)*(z9-zcm) + (z10-zcm)*(z10-zcm))/10 + +# Calculating the eigen values of the gyration tensor +min_eig = [] + +for i in range(len(t)): + G = np.matrix([ [G11[i],G12[i],G13[i]], + [G21[i],G22[i],G23[i]], + [G31[i],G32[i],G33[i]] ]) + + eigenvalues = np.linalg.eigvals(G) + eig = np.min(eigenvalues) + min_eig.append(eig) + +plt.plot(t,min_eig,linewidth = 1) +plt.xlim([0,0.5]) +plt.xlabel('Time (s)') +plt.ylabel('Min. eigen value of G') +plt.minorticks_on() +plt.grid() +plt.savefig('RBC_Posieuille_Gyration.eps') +plt.show() + + +t = t*1000 +yr = ycm/10 + +plt.plot(zcm,ycm) +plt.show() + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh new file mode 100755 index 00000000..7bd597b8 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh @@ -0,0 +1,31 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allclean script for redBloodCellShearFlow test case +# clean up redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#------------------------------------------------------------------------------ +#- clean up case +echo "deleting data at: $casePath" +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanTimeDirectories +rm -rf processor* > /dev/null 2>&1 +rm -rf clockData > /dev/null 2>&1 +echo "Cleaning $casePath/DEM/post" +rm $casePath/DEM/post/*.* > /dev/null 2>&1 +#rm -r $casePath/DEM/post/restart/*.* +#touch $casePath/DEM/post/restart/.gitignore +echo "done" + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh new file mode 100755 index 00000000..00761fb1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh @@ -0,0 +1,30 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allrun script for redBloodCellShearFlow test case +# run redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +# check if mesh was built +if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then + echo "mesh was built before - using old mesh" +else + echo "mesh needs to be built" + cd $casePath/CFD + blockMesh +fi + +if [ -f "$casePath/DEM/post/restart/liggghts.restart" ]; then + echo "LIGGGHTS init was run before - using existing restart file" +else + #- run serial DEM + $casePath/DEMrun.sh +fi + +#- run parallel CFD-DEM +bash $casePath/parCFDDEMrun.sh diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U new file mode 100644 index 00000000..4c30855a --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); //velocity in micrometres per second + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type fixedValue; + value uniform ( 0.00225 0 0); + } + bottom + { + type fixedValue; + value uniform (-0.00225 0 0); + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us new file mode 100644 index 00000000..9fadb65e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object Us; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda new file mode 100644 index 00000000..8d3e3b8f --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object lambda; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 100; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p new file mode 100644 index 00000000..c18c6019 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB new file mode 100644 index 00000000..30eb5d9e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object phiIB; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction new file mode 100644 index 00000000..e1bb7586 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object voidfraction; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties new file mode 100644 index 00000000..b156d097 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties @@ -0,0 +1,92 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object couplingProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// sub-models & settings + +modelType none; + +checkPeriodicCells; + +couplingInterval 10; + +depth 0; + +voidFractionModel IB; + +locateModel engineIB; + +meshMotionModel noMeshMotion; + +regionModel allRegion; + +dataExchangeModel twoWayMPI; + +IOModel off; + +probeModel off; + +averagingModel dilute; + +clockModel off; + +smoothingModel off; + +forceModels +( + ShirgaonkarIB +); + +momCoupleModels +( +); + +turbulenceModelType "turbulenceProperties"; + + +// sub-model properties + +ShirgaonkarIBProps +{ + velFieldName "U"; + pressureFieldName "p"; + densityFieldName "rho"; + solidVolFractionName "lambda"; + useTorque; +} + +IBProps +{ + maxCellsPerParticle 20000; + alphaMin 0.30; + scaleUpVol 1.0; +} + +engineIBProps +{ + engineProps + { + treeSearch true; + } + zSplit 8; + xySplit 16; +} + +twoWayMPIProps +{ + liggghtsPath "../DEM/in.liggghts_run"; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g new file mode 100644 index 00000000..d4b9882c --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 0 0); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands new file mode 100644 index 00000000..9a3b6ab0 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object liggghtsCommands; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +liggghtsCommandModels +( + runLiggghts +); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties new file mode 100644 index 00000000..80258c87 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +transportModel Newtonian; + +// Fluid properties of plasma // dextran +nu nu [ 0 2 -1 0 0 0 0 ] 1.2; //4.9524; + +CrossPowerLawCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + m m [ 0 0 1 0 0 0 0 ] 1; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +BirdCarreauCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [ 0 0 1 0 0 0 0 ] 0; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties new file mode 100644 index 00000000..f6d9d4cb --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict new file mode 100644 index 00000000..e34cf832 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict @@ -0,0 +1,110 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + ( 0.0 0.0 0.0) + (45.0 0.0 0.0) + (45.0 22.5 0.0) + ( 0.0 22.5 0.0) + ( 0.0 0.0 22.5) + (45.0 0.0 22.5) + (45.0 22.5 22.5) + ( 0.0 22.5 22.5) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (58 29 29) simpleGrading (1 1 1) +); + +boundary +( + front + { + type cyclicAMI; + neighbourPatch back; + faces + ( + (4 5 6 7) + ); + transform translational; + separationVector (0 0 -22.5); + } + + back + { + type cyclicAMI; + neighbourPatch front; + faces + ( + (3 2 1 0) + ); + transform translational; + separationVector (0 0 22.5); + } + + top + { + type wall; + faces + ( + (2 3 7 6) + ); + } + + bottom + { + type wall; + faces + ( + (0 1 5 4) + ); + } + + left + { + type cyclicAMI; + neighbourPatch right; + faces + ( + (0 4 7 3) + ); + transform translational; + separationVector ( 45 0 0); + } + + right + { + type cyclicAMI; + neighbourPatch left; + faces + ( + (1 2 6 5) + ); + transform translational; + separationVector (-45 0 0); + } +); + +mergePatchPairs +( +); + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict new file mode 100644 index 00000000..f86ad4d8 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict @@ -0,0 +1,121 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application pisoFoam; + +startFrom latestTime; + +startTime 0; + +stopAt endTime; + +endTime 100000; + +deltaT 1; + +writeControl adjustableRunTime; + +writeInterval 100; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep no; + +maxCo 1.0; + +maxDeltaT 10; + + +// ************************************************************************* // +DimensionedConstants +{ + unitSet micro; // SI; + + SICoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-11; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-34; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-31; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-27; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-27; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } + + microCoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-20; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-13; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-16; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-12; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-12; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } +} diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict new file mode 100644 index 00000000..6f0e5969 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict @@ -0,0 +1,29 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 2; + +method simple; + +simpleCoeffs +{ + n (1 1 2); + delta 0.001; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions new file mode 100644 index 00000000..e41aadb1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object fvOptions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +momentumSource +{ + type meanVelocityForce; + active no; + + meanVelocityForceCoeffs + { + selectionMode all; + fields (U); + // Note: set mean velocity here + Ubar ( 0 0 0.0004 ); + relaxation 0.9; + } +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes new file mode 100644 index 00000000..005e8180 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss skewCorrected linear; + grad(p) Gauss skewCorrected linear; + grad(U) Gauss skewCorrected linear; +} + +divSchemes +{ + default Gauss skewCorrected linear; + div(phi,U) Gauss skewCorrected upwind; + div(phi,k) Gauss limitedLinear 1; + div(phi,epsilon) Gauss limitedLinear 1; + div(phi,R) Gauss limitedLinear 1; + div(R) Gauss linear; + div(phi,nuTilda) Gauss limitedLinear 1; + div((nuEff*dev(grad(U).T()))) Gauss linear; + div(U) Gauss skewCorrected lineaUpwind; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(nuEff,U) Gauss linear corrected; + laplacian((1|A(U)),p) Gauss linear corrected; + laplacian((voidfraction2|A(U)),p) Gauss linear corrected; + laplacian(DkEff,k) Gauss linear corrected; + laplacian(DepsilonEff,epsilon) Gauss linear corrected; + laplacian(DREff,R) Gauss linear corrected; + laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; + laplacian(phiIB) Gauss linear corrected; + laplacian(U) Gauss linear corrected; +} + +interpolationSchemes +{ + default skewCorrected linear; + interpolate(U) skewCorrected linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution new file mode 100644 index 00000000..b4c9e0ec --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution @@ -0,0 +1,118 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + pFinal + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + U + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + 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; + } + + phiIB + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + "(U|k)Final" + { + $U; + tolerance 1e-05; + relTol 0; + } + + Us + { + $U; + tolerance 1e-05; + relTol 0; + } + +} + +PISO +{ + nCorrectors 2; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + +relaxationFactors +{ +/* + equations + { + ".*" 0.9; + } +*/ +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt new file mode 100644 index 00000000..55c1fc96 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt @@ -0,0 +1,10 @@ +1.58162822 1.14912017 0.00000000 1.95500000 +0.60412822 1.85931549 0.00000000 1.95500000 +-0.60412822 1.85931549 0.00000000 1.95500000 +-1.58162822 1.14912017 0.00000000 1.95500000 +-1.95500000 0.00000000 0.00000000 1.95500000 +-1.58162822 -1.14912017 0.00000000 1.95500000 +-0.60412822 -1.85931549 0.00000000 1.95500000 +0.60412822 -1.85931549 0.00000000 1.95500000 +1.58162822 -1.14912017 0.00000000 1.95500000 +1.95500000 -0.00000000 0.00000000 1.95500000 diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init new file mode 100644 index 00000000..c1fb1d62 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init @@ -0,0 +1,113 @@ +# Insert a RBC into a box, then induce shear flow +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 300 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Sben equal ${Sn}*10 # analogous bending modulus of the parallel bond +variable Stor equal ${Sn}*10 # analogous torsional modulus of the parallel bond + + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary p f p + +# communication +newton off +communicate single vel yes + +# units +units micro + +#region and box creation +region domain block 0 45 0 22.5 0 22.5 units box +create_box 2 domain + +# neighbors +neighbor 1 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e03 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.5 +fix m3 all property/global coefficientRestitution peratomtypepair 2 1e-06 1e-03 1e-03 1e-06 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.1 0.1 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 2.65 0.0 1 1e50 1e50 + + +# particle template, insertion region and bond formation +fix pts1 all particletemplate/multiplespheres 15485863 atom_type 1 & +density constant 1098 nspheres 10 ntry 1000000 spheres file ../DEM/data/spheres.txt scale 1 & +bonded yes/explicit & +nbond_pairs 10 & +1 2 & +2 3 & +3 4 & +4 5 & +5 6 & +6 7 & +7 8 & +8 9 & +9 10 & +10 1 & +bond_type 1 + +fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0 + +region bc1 block 0 45 0 22.5 0 22.5 units box + +fix ins1 all insert/pack seed 32452867 distributiontemplate pdd1 & + maxattempt 10000 insert_every once overlapcheck yes vel constant 0.0 0.0 0.0 & + region bc1 particles_in_region 1 pos 22.5 11.25 11.25 ntry_mc 100000 + +run 1000 + +# set timestep +timestep 0.001 + +# check timestep +fix ts_check all check/timestep/gran 1000 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 10000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump*.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bondtype f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds*.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 100000 + +write_restart ../DEM/post/restart/liggghts.restart + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run new file mode 100644 index 00000000..822edd9e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run @@ -0,0 +1,130 @@ +# RBC in a box sheared by fluid +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 25 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Ab equal (22/7)*${rb}${rb} # bond area +variable Sben equal ${Sn}/${Ab} # analogous bending modulus of the parallel bond +variable Stor equal ${St}/${Ab} # analogous torsional modulus of the parallel bond + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary p f p + +# communication +newton off +communicate single vel yes + +# units and processors +units micro +processors 1 1 2 + +# read the restart file +read_restart ../DEM/post/restart/liggghts.restart + +# neighbors +neighbor 1.0 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e04 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.18 +fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.1 0.1 0.9 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.3 0.3 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 3.0 1.0 1.0 1e50 1e50 + + +# reset timestep after reading the restart file +reset_timestep 0 +timestep 0.1 + +# cfd coupling +fix cfd all couple/cfd couple_every 10 mpi +fix cfd2 all couple/cfd/force + +# apply nve integration to all constituent spheres - single sphere treatment +fix integr all nve + +compute unwra all property/atom xu yu zu + +# variables for data processing +variable x1 equal c_unwra[1][1] +variable y1 equal c_unwra[1][2] +variable z1 equal c_unwra[1][3] +variable x2 equal c_unwra[2][1] +variable y2 equal c_unwra[2][2] +variable z2 equal c_unwra[2][3] +variable x3 equal c_unwra[3][1] +variable y3 equal c_unwra[3][2] +variable z3 equal c_unwra[3][3] +variable x4 equal c_unwra[4][1] +variable y4 equal c_unwra[4][2] +variable z4 equal c_unwra[4][3] +variable x5 equal c_unwra[5][1] +variable y5 equal c_unwra[5][2] +variable z5 equal c_unwra[5][3] +variable x6 equal c_unwra[6][1] +variable y6 equal c_unwra[6][2] +variable z6 equal c_unwra[6][3] +variable x7 equal c_unwra[7][1] +variable y7 equal c_unwra[7][2] +variable z7 equal c_unwra[7][3] +variable x8 equal c_unwra[8][1] +variable y8 equal c_unwra[8][2] +variable z8 equal c_unwra[8][3] +variable x9 equal c_unwra[9][1] +variable y9 equal c_unwra[9][2] +variable z9 equal c_unwra[9][3] +variable x10 equal c_unwra[10][1] +variable y10 equal c_unwra[10][2] +variable z10 equal c_unwra[10][3] +variable time equal step*dt + +fix extra1 all print 10 "${time} ${x1} ${y1} ${z1} ${x2} ${y2} ${z2} ${x3} ${y3} ${z3} ${x4} ${y4} ${z4} ${x5} ${y5} ${z5} ${x6} ${y6} ${z6} ${x7} ${y7} ${z7} ${x8} ${y8} ${z8} ${x9} ${y9} ${z9} ${x10} ${y10} ${z10}" & + file ../DEM/post/particle_position.txt screen no + +# check timestep +fix ts_check all check/timestep/gran 1 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 1000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bondtype f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 1 + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/post/restart/.gitignore b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/post/restart/.gitignore new file mode 100644 index 00000000..e69de29b diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh new file mode 100755 index 00000000..94bd244f --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh @@ -0,0 +1,26 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# DEMrun script for redBloodCellShearFlow test case +# init redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- include functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath="$casePath" +headerText="run_liggghts_init_DEM" +logfileName="log_$headerText" +solverName="in.liggghts_init" +debugMode="off" +#------------------------------------------------------------------------------ + +#- call function to run DEM case +DEMrun $logpath $logfileName $casePath $headerText $solverName $debugMode + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md new file mode 100644 index 00000000..94373848 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md @@ -0,0 +1,9 @@ +## Single Red Blood Cell Dynamics in Shear Flow + +This case uses a continuous forcing immersed boundary solver on the CFD side and +a reduced order model of a red blood cell on the DEM side, cf. + +Balachandran Nair, A.N., Pirker, S. and Saeedipour, M., "Resolved CFD-DEM +simulation of blood flow with a reduced-order RBC model", Comp. Part. Mech. +(2021) + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh new file mode 100755 index 00000000..a0cd6f60 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh @@ -0,0 +1,59 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# parCFDDEMrun script for redBloodCellShearFlow test case +# run redBloodCellShearFlow CFD-DEM +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath=$casePath +headerText="run_parallel_cfdemSolverIBContinuousForcing_redBloodCellShearFlow" +logfileName="log_$headerText" +solverName="cfdemSolverIBContinuousForcing" +nrProcs="2" +machineFileName="none" # yourMachinefileName | none +debugMode="off" # on | off| strict +testHarnessPath="$CFDEM_TEST_HARNESS_PATH" +runPython="false" +postproc="false" +#------------------------------------------------------------------------------ + +#- call function to run a parallel CFD-DEM case +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode + +if [ $runPython == "true" ] + then + cd $casePath + python results.py +fi + +if [ $postproc == "true" ] + then + #- get VTK data from liggghts dump file + cd $casePath/DEM/post + lpp dump* + + #- get VTK data from CFD sim + cd $casePath/CFD + reconstructParMesh -constant -mergeTol 1e-06 + reconstructPar -noLagrangian + foamToVTK + + #- start paraview + paraview + + #- keep terminal open (if started in new terminal) + echo "...press enter to clean up case" + echo "press Ctr+C to keep data" + read +fi + + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh new file mode 100755 index 00000000..a3efb45b --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh @@ -0,0 +1,31 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# post script for redBloodCellShearFlow test case +# run post-processing +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +source ~/.bashrc + +#- source CFDEM functions +shopt -s expand_aliases +source $CFDEM_PROJECT_DIR/etc/bashrc + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#------------------------------------------------------------------------------ +#- get VTK data from liggghts dump file +cd $casePath/DEM/post +lpp dump* + +#- get VTK data from CFD sim +cd $casePath/CFD +reconstructParMesh -constant -mergeTol 1e-06 +reconstructPar -noLagrangian + +#- start paraview +paraview + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py new file mode 100755 index 00000000..f11fbfc1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py @@ -0,0 +1,66 @@ +#!/usr/bin/python + +import numpy as np +import matplotlib.pyplot as plt + +# Reading sphere data file +t,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,x9,y9,z9,x10,y10,z10 = np.genfromtxt('DEM/post/particle_position.txt', skip_header=0,usecols = (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30),unpack=True) +t = t*1e-06 + +# Calculating the centre of mass of RBC +xcm = (x1+x2+x3+x4+x5+x6+x7+x8+x9+x10)/10 +ycm = (y1+y2+y3+y4+y5+y6+y7+y8+y9+y10)/10 +zcm = (z1+z2+z3+z4+z5+z6+z7+z8+z9+z10)/10 + +# Calculating the gyration tensor +G11 = ((x1-xcm)*(x1-xcm) + (x2-xcm)*(x2-xcm) + (x3-xcm)*(x3-xcm) + (x4-xcm)*(x4-xcm) + (x5-xcm)*(x5-xcm) + (x6-xcm)*(x6-xcm) + (x7-xcm)*(x7-xcm) + (x8-xcm)*(x8-xcm) + (x9-xcm)*(x9-xcm) + (x10-xcm)*(x10-xcm))/10 +G12 = ((x1-xcm)*(y1-ycm) + (x2-xcm)*(y2-ycm) + (x3-xcm)*(y3-ycm) + (x4-xcm)*(y4-ycm) + (x5-xcm)*(y5-ycm) + (x6-xcm)*(y6-ycm) + (x7-xcm)*(y7-ycm) + (x8-xcm)*(y8-ycm) + (x9-xcm)*(y9-ycm) + (x10-xcm)*(y10-ycm))/10 +G13 = ((x1-xcm)*(z1-zcm) + (x2-xcm)*(z2-zcm) + (x3-xcm)*(z3-zcm) + (x4-xcm)*(z4-zcm) + (x5-xcm)*(z5-zcm) + (x6-xcm)*(z6-zcm) + (x7-xcm)*(z7-zcm) + (x8-xcm)*(z8-zcm) + (x9-xcm)*(z9-zcm) + (x10-xcm)*(z10-zcm))/10 +G21 = ((y1-ycm)*(x1-xcm) + (y2-ycm)*(x2-xcm) + (y3-ycm)*(x3-xcm) + (y4-ycm)*(x4-xcm) + (y5-ycm)*(x5-xcm) + (y6-ycm)*(x6-xcm) + (y7-ycm)*(x7-xcm) + (y8-ycm)*(x8-xcm) + (y9-ycm)*(x9-xcm) + (y10-ycm)*(x10-xcm))/10 +G22 = ((y1-ycm)*(y1-ycm) + (y2-ycm)*(y2-ycm) + (y3-ycm)*(y3-ycm) + (y4-ycm)*(y4-ycm) + (y5-ycm)*(y5-ycm) + (y6-ycm)*(y6-ycm) + (y7-ycm)*(y7-ycm) + (y8-ycm)*(y8-ycm) + (y9-ycm)*(y9-ycm) + (y10-ycm)*(y10-ycm))/10 +G23 = ((y1-ycm)*(z1-zcm) + (y2-ycm)*(z2-zcm) + (y3-ycm)*(z3-zcm) + (y4-ycm)*(z4-zcm) + (y5-ycm)*(z5-zcm) + (y6-ycm)*(z6-zcm) + (y7-ycm)*(z7-zcm) + (y8-ycm)*(z8-zcm) + (y9-ycm)*(z9-zcm) + (y10-ycm)*(z10-zcm))/10 +G31 = ((z1-zcm)*(x1-xcm) + (z2-zcm)*(x2-xcm) + (z3-zcm)*(x3-xcm) + (z4-zcm)*(x4-xcm) + (z5-zcm)*(x5-xcm) + (z6-zcm)*(x6-xcm) + (z7-zcm)*(x7-xcm) + (z8-zcm)*(x8-xcm) + (z9-zcm)*(x9-xcm) + (z10-zcm)*(x10-xcm))/10 +G32 = ((z1-zcm)*(y1-ycm) + (z2-zcm)*(y2-ycm) + (z3-zcm)*(y3-ycm) + (z4-zcm)*(y4-ycm) + (z5-zcm)*(y5-ycm) + (z6-zcm)*(y6-ycm) + (z7-zcm)*(y7-ycm) + (z8-zcm)*(y8-ycm) + (z9-zcm)*(y9-ycm) + (z10-zcm)*(y10-ycm))/10 +G33 = ((z1-zcm)*(z1-zcm) + (z2-zcm)*(z2-zcm) + (z3-zcm)*(z3-zcm) + (z4-zcm)*(z4-zcm) + (z5-zcm)*(z5-zcm) + (z6-zcm)*(z6-zcm) + (z7-zcm)*(z7-zcm) + (z8-zcm)*(z8-zcm) + (z9-zcm)*(z9-zcm) + (z10-zcm)*(z10-zcm))/10 + +# Calculating the eigen values of the gyration tensor +min_eig = [] + +for i in range(len(t)): + G = np.matrix([ [G11[i],G12[i],G13[i]], + [G21[i],G22[i],G23[i]], + [G31[i],G32[i],G33[i]] ]) + + eigenvalues = np.linalg.eigvals(G) + eig = np.min(eigenvalues) + min_eig.append(eig) + +plt.plot(t,min_eig,linewidth = 2) +plt.xlabel('Time (s)') +plt.ylabel('Min. eigen value of G') +plt.grid() +plt.show() + +plt.plot(z1,y1) +plt.plot(z2,y2) +plt.plot(z3,y3) +plt.plot(z4,y4) +plt.plot(z5,y5) +plt.plot(z6,y6) +plt.plot(z7,y7) +plt.plot(z8,y8) +plt.plot(z9,y9) +plt.plot(z10,y10) +plt.show() + +plt.plot(z1,x1) +plt.plot(z2,x2) +plt.plot(z3,x3) +plt.plot(z4,x4) +plt.plot(z5,x5) +plt.plot(z6,x6) +plt.plot(z7,x7) +plt.plot(z8,x8) +plt.plot(z9,x9) +plt.plot(z10,x10) +plt.show()