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

Merge pull request #130 from ParticulateFlow/feature/immersed_multisphere

Browse Source 
close #127
This commit is contained in:
Daniel
2022-03-25 09:14:06 +01:00
committed by GitHub
parent 7698e20c8f 8cb7184a7c
commit 32e41e2fa6
97 changed files with 4635 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1
doc/CFDEMcoupling_solvers.txt
View File

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

82
doc/cfdemSolverIBContinuousForcing.txt Normal file
View File

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

2
doc/forceModel_ShirgaonkarIB.txt
View File

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

1
etc/solver-list.txt
View File

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

2
etc/tutorial-list.txt
View File

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

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

@ -36,7 +36,6 @@ Description
#include "locateModel.H"
#include "dataExchangeModel.H"
#include "IOModel.H"
#include <mpi.h>
#include "IOmanip.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -58,6 +57,7 @@ cfdemCloudIB::cfdemCloudIB
angularVelocities_(NULL),
pRefCell_(readLabel(mesh.solutionDict().subDict("PISO").lookup("pRefCell"))),
pRefValue_(readScalar(mesh.solutionDict().subDict("PISO").lookup("pRefValue"))),
DEMTorques_(NULL),
haveEvolvedOnce_(false),
skipLagrangeToEulerMapping_(false)
{
@ -75,6 +75,7 @@ cfdemCloudIB::cfdemCloudIB
cfdemCloudIB::~cfdemCloudIB()
{
dataExchangeM().destroy(angularVelocities_,3);
dataExchangeM().destroy(DEMTorques_,3);
}
@ -82,8 +83,7 @@ cfdemCloudIB::~cfdemCloudIB()
void cfdemCloudIB::getDEMdata()
{
cfdemCloud::getDEMdata();
Info << "=== cfdemCloudIB::getDEMdata() === particle rotation not considered in CFD" << endl;
//dataExchangeM().getData("omega","vector-atom",angularVelocities_);
dataExchangeM().getData("omega","vector-atom",angularVelocities_);
}
bool cfdemCloudIB::reAllocArrays()
@ -92,12 +92,24 @@ bool cfdemCloudIB::reAllocArrays()
{
// get arrays of new length
dataExchangeM().allocateArray(angularVelocities_,0.,3);
dataExchangeM().allocateArray(DEMTorques_,0.,3);
return true;
}
return false;
}
bool cfdemCloudIB::evolve()
void cfdemCloudIB::giveDEMdata()
{
cfdemCloud::giveDEMdata();
dataExchangeM().giveData("hdtorque","vector-atom",DEMTorques_);
}
inline double ** cfdemCloudIB::DEMTorques() const
{
return DEMTorques_;
}
bool cfdemCloudIB::evolve(volVectorField& Us)
{
numberOfParticlesChanged_ = false;
arraysReallocated_=false;
@ -140,6 +152,10 @@ bool cfdemCloudIB::evolve()
for (int i=0;i<nrForceModels();i++) forceM(i).setForce();
if(verbose_) Info << "setForce done." << endl;
// set particle velocity field
if(verbose_) Info << "- setVelocity(velocities_)" << endl;
calcForcingTerm(Us); // only needed for continuous forcing
// write DEM data
if(verbose_) Info << " -giveDEMdata()" << endl;
giveDEMdata();
@ -158,6 +174,34 @@ bool cfdemCloudIB::evolve()
return doCouple;
}
void cfdemCloudIB::calcForcingTerm(volVectorField& Us)
{
label cell = 0;
vector uP(0,0,0);
vector rVec(0,0,0);
vector velRot(0,0,0);
vector angVel(0,0,0);
Us.primitiveFieldRef() = vector::zero;
for (int index = 0; index < numberOfParticles(); ++index)
{
for (int subCell = 0; subCell < voidFractionM().cellsPerParticle()[index][0]; subCell++)
{
cell = cellIDs()[index][subCell];
if (cell >=0)
{
// calc particle velocity
for (int i=0;i<3;i++) rVec[i]=Us.mesh().C()[cell][i]-position(index)[i];
for (int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i];
velRot=angVel^rVec;
for (int i=0;i<3;i++) uP[i] = velocities()[index][i]+velRot[i];
Us[cell] = (1-voidfractions_[index][subCell])*uP;
}
}
}
Us.correctBoundaryConditions();
}
void cfdemCloudIB::calcVelocityCorrection
(
volScalarField& p,
@ -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"))

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@ -0,0 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class 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;
}
}
// ************************************************************************* //

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

@ -0,0 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class 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;
}
}
// ************************************************************************* //

39
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho Normal file
View File

@ -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;
}
}
// ************************************************************************* //

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

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

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

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

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

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

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

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

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

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

40
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/transportProperties Executable file
View File

@ -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;
}
// ************************************************************************* //

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

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

33
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/octave/plot_data.m Normal file
View File

@ -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')

77
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/blockMeshDict Normal file
View File

@ -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
(
);
// ************************************************************************* //

56
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/controlDict Normal file
View File

@ -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;
// ************************************************************************* //

32
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/decomposeParDict Normal file
View File

@ -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;
}
// ************************************************************************* //

34
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvOptions Normal file
View File

@ -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;
}
}
// ************************************************************************* //

75
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSchemes Executable file
View File

@ -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 ;
}
// ************************************************************************* //

110
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/system/fvSolution Executable file
View File

@ -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;
}
// ************************************************************************* //

88
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/in.liggghts_run Normal file
View File

@ -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

0
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/DEM/post/.gitignore vendored Normal file
View File

74
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/parCFDDEMrun.sh Normal file
View File

@ -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"

32
tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/post.sh Executable file
View File

@ -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

31
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allclean.sh Executable file
View File

@ -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"

30
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/Allrun.sh Executable file
View File

@ -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

39
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/U Normal file
View File

@ -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;
}
}
// ************************************************************************* //

42
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/Us Normal file
View File

@ -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;
}
}
// ************************************************************************* //

40
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/lambda Normal file
View File

@ -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;
}
}
// ************************************************************************* //

40
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/p Normal file
View File

@ -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;
}
}
// ************************************************************************* //

41
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/phiIB Normal file
View File

@ -0,0 +1,41 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class 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;
}
}
// ************************************************************************* //

39
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/0/voidfraction Normal file
View File

@ -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;
}
}
// ************************************************************************* //

91
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/couplingProperties Normal file
View File

@ -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";
}
// ************************************************************************* //

21
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/g Normal file
View File

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

22
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/liggghtsCommands Normal file
View File

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

39
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/transportProperties Normal file
View File

@ -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;
}
// ************************************************************************* //

20
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/constant/turbulenceProperties Normal file
View File

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

123
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/blockMeshDict Normal file
View File

@ -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
(
);
// ************************************************************************* //

121
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/controlDict Normal file
View File

@ -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;
}
}
}

28
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/decomposeParDict Normal file
View File

@ -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;
}
// ************************************************************************* //

34
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvOptions Normal file
View File

@ -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;
}
}
// ************************************************************************* //

75
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSchemes Normal file
View File

@ -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 ;
}
// ************************************************************************* //

116
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/CFD/system/fvSolution Normal file
View File

@ -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;
}
}
// ************************************************************************* //

10
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/data/spheres.txt Normal file
View File

@ -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

116
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_init Normal file
View File

@ -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

133
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/in.liggghts_run Normal file
View File

@ -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

0
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEM/post/restart/.gitignore vendored Normal file
View File

26
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/DEMrun.sh Executable file
View File

@ -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

59
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/parCFDDEMrun.sh Executable file
View File

@ -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

57
tutorials/cfdemSolverIBContinuousForcing/redBloodCellPoiseuilleFlow/results.py Executable file
View File

@ -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()

31
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh Executable file
View File

@ -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"

30
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh Executable file
View File

@ -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

51
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U Normal file
View File

@ -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;
}
}
// ************************************************************************* //

50
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us Normal file
View File

@ -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;
}
}
// ************************************************************************* //

49
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda Normal file
View File

@ -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;
}
}
// ************************************************************************* //

49
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p Normal file
View File

@ -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;
}
}
// ************************************************************************* //

49
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB Normal file
View File

@ -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;
}
}
// ************************************************************************* //

49
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction Normal file
View File

@ -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;
}
}
// ************************************************************************* //

92
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties Normal file
View File

@ -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";
}
// ************************************************************************* //

21
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g Normal file
View File

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

22
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands Normal file
View File

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

40
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties Normal file
View File

@ -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;
}
// ************************************************************************* //

20
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties Normal file
View File

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

110
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict Normal file
View File

@ -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
(
);
// ************************************************************************* //

121
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict Normal file
View File

@ -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;
}
}
}

29
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict Normal file
View File

@ -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;
}
// ************************************************************************* //

34
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions Normal file
View File

@ -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;
}
}
// ************************************************************************* //

75
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes Normal file
View File

@ -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 ;
}
// ************************************************************************* //

118
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution Normal file
View File

@ -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;
}
*/
}
// ************************************************************************* //

10
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt Normal file
View File

@ -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

113
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init Normal file
View File

@ -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

130
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run Normal file
View File

@ -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

0
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/post/restart/.gitignore vendored Normal file
View File

26
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh Executable file
View File

@ -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

9
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md Normal file
View File

@ -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)

59
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh Executable file
View File

@ -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

31
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh Executable file
View File

@ -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

66
tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py Executable file
View File

@ -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()