add a tutorial for cfdemSolverIBContinuousForcing

simulating the deformation of a red blood cell (bonded particles) in
shear flow
This commit is contained in:
danielque
2022-03-24 15:39:10 +01:00
parent 38b8d6c8b8
commit 1f5c8f6492
28 changed files with 1484 additions and 0 deletions

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"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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)

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

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

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