diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh new file mode 100755 index 00000000..7bd597b8 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allclean.sh @@ -0,0 +1,31 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allclean script for redBloodCellShearFlow test case +# clean up redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#------------------------------------------------------------------------------ +#- clean up case +echo "deleting data at: $casePath" +source $WM_PROJECT_DIR/bin/tools/CleanFunctions +cd $casePath/CFD +cleanTimeDirectories +rm -rf processor* > /dev/null 2>&1 +rm -rf clockData > /dev/null 2>&1 +echo "Cleaning $casePath/DEM/post" +rm $casePath/DEM/post/*.* > /dev/null 2>&1 +#rm -r $casePath/DEM/post/restart/*.* +#touch $casePath/DEM/post/restart/.gitignore +echo "done" + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh new file mode 100755 index 00000000..00761fb1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/Allrun.sh @@ -0,0 +1,30 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# Allrun script for redBloodCellShearFlow test case +# run redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +# check if mesh was built +if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then + echo "mesh was built before - using old mesh" +else + echo "mesh needs to be built" + cd $casePath/CFD + blockMesh +fi + +if [ -f "$casePath/DEM/post/restart/liggghts.restart" ]; then + echo "LIGGGHTS init was run before - using existing restart file" +else + #- run serial DEM + $casePath/DEMrun.sh +fi + +#- run parallel CFD-DEM +bash $casePath/parCFDDEMrun.sh diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U new file mode 100644 index 00000000..4c30855a --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/U @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); //velocity in micrometres per second + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type fixedValue; + value uniform ( 0.00225 0 0); + } + bottom + { + type fixedValue; + value uniform (-0.00225 0 0); + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us new file mode 100644 index 00000000..9fadb65e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/Us @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object Us; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda new file mode 100644 index 00000000..8d3e3b8f --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/lambda @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object lambda; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 100; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p new file mode 100644 index 00000000..c18c6019 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/p @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB new file mode 100644 index 00000000..30eb5d9e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/phiIB @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object phiIB; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction new file mode 100644 index 00000000..e1bb7586 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/0/voidfraction @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object voidfraction; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + front + { + type cyclicAMI; + } + back + { + type cyclicAMI; + } + top + { + type zeroGradient; + } + bottom + { + type zeroGradient; + } + left + { + type cyclicAMI; + } + right + { + type cyclicAMI; + } +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties new file mode 100644 index 00000000..b156d097 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/couplingProperties @@ -0,0 +1,92 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object couplingProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// sub-models & settings + +modelType none; + +checkPeriodicCells; + +couplingInterval 10; + +depth 0; + +voidFractionModel IB; + +locateModel engineIB; + +meshMotionModel noMeshMotion; + +regionModel allRegion; + +dataExchangeModel twoWayMPI; + +IOModel off; + +probeModel off; + +averagingModel dilute; + +clockModel off; + +smoothingModel off; + +forceModels +( + ShirgaonkarIB +); + +momCoupleModels +( +); + +turbulenceModelType "turbulenceProperties"; + + +// sub-model properties + +ShirgaonkarIBProps +{ + velFieldName "U"; + pressureFieldName "p"; + densityFieldName "rho"; + solidVolFractionName "lambda"; + useTorque; +} + +IBProps +{ + maxCellsPerParticle 20000; + alphaMin 0.30; + scaleUpVol 1.0; +} + +engineIBProps +{ + engineProps + { + treeSearch true; + } + zSplit 8; + xySplit 16; +} + +twoWayMPIProps +{ + liggghtsPath "../DEM/in.liggghts_run"; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g new file mode 100644 index 00000000..d4b9882c --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 0 0); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands new file mode 100644 index 00000000..9a3b6ab0 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/liggghtsCommands @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object liggghtsCommands; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +liggghtsCommandModels +( + runLiggghts +); + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties new file mode 100644 index 00000000..80258c87 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/transportProperties @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +transportModel Newtonian; + +// Fluid properties of plasma // dextran +nu nu [ 0 2 -1 0 0 0 0 ] 1.2; //4.9524; + +CrossPowerLawCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + m m [ 0 0 1 0 0 0 0 ] 1; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +BirdCarreauCoeffs +{ + nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06; + nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [ 0 0 1 0 0 0 0 ] 0; + n n [ 0 0 0 0 0 0 0 ] 1; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties new file mode 100644 index 00000000..f6d9d4cb --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/constant/turbulenceProperties @@ -0,0 +1,20 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict new file mode 100644 index 00000000..e34cf832 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/blockMeshDict @@ -0,0 +1,110 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + ( 0.0 0.0 0.0) + (45.0 0.0 0.0) + (45.0 22.5 0.0) + ( 0.0 22.5 0.0) + ( 0.0 0.0 22.5) + (45.0 0.0 22.5) + (45.0 22.5 22.5) + ( 0.0 22.5 22.5) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (58 29 29) simpleGrading (1 1 1) +); + +boundary +( + front + { + type cyclicAMI; + neighbourPatch back; + faces + ( + (4 5 6 7) + ); + transform translational; + separationVector (0 0 -22.5); + } + + back + { + type cyclicAMI; + neighbourPatch front; + faces + ( + (3 2 1 0) + ); + transform translational; + separationVector (0 0 22.5); + } + + top + { + type wall; + faces + ( + (2 3 7 6) + ); + } + + bottom + { + type wall; + faces + ( + (0 1 5 4) + ); + } + + left + { + type cyclicAMI; + neighbourPatch right; + faces + ( + (0 4 7 3) + ); + transform translational; + separationVector ( 45 0 0); + } + + right + { + type cyclicAMI; + neighbourPatch left; + faces + ( + (1 2 6 5) + ); + transform translational; + separationVector (-45 0 0); + } +); + +mergePatchPairs +( +); + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict new file mode 100644 index 00000000..f86ad4d8 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/controlDict @@ -0,0 +1,121 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application pisoFoam; + +startFrom latestTime; + +startTime 0; + +stopAt endTime; + +endTime 100000; + +deltaT 1; + +writeControl adjustableRunTime; + +writeInterval 100; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep no; + +maxCo 1.0; + +maxDeltaT 10; + + +// ************************************************************************* // +DimensionedConstants +{ + unitSet micro; // SI; + + SICoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-11; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-34; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-31; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-27; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-27; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } + + microCoeffs + { + universal + { + c c [ 0 1 -1 0 0 0 0 ] 2.99792e+08; + G G [ -1 3 -2 0 0 0 0 ] 6.67429e-20; + h h [ 1 2 -1 0 0 0 0 ] 6.62607e-13; + } + electromagnetic + { + e e [ 0 0 1 0 0 1 0 ] 1.60218e-19; + } + atomic + { + me me [ 1 0 0 0 0 0 0 ] 9.10938e-16; + mp mp [ 1 0 0 0 0 0 0 ] 1.67262e-12; + } + physicoChemical + { + mu mu [ 1 0 0 0 0 0 0 ] 1.66054e-12; + k k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23; + } + standard + { + //- Standard pressure [Pa] + Pstd Pstd [ 1 -1 -2 0 0 0 0 ] 100000; + //- Standard temperature [degK] + Tstd Tstd [ 0 0 0 1 0 0 0 ] 298.15; + } + } +} diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict new file mode 100644 index 00000000..6f0e5969 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/decomposeParDict @@ -0,0 +1,29 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 2; + +method simple; + +simpleCoeffs +{ + n (1 1 2); + delta 0.001; +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions new file mode 100644 index 00000000..e41aadb1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvOptions @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object fvOptions; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +momentumSource +{ + type meanVelocityForce; + active no; + + meanVelocityForceCoeffs + { + selectionMode all; + fields (U); + // Note: set mean velocity here + Ubar ( 0 0 0.0004 ); + relaxation 0.9; + } +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes new file mode 100644 index 00000000..005e8180 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSchemes @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss skewCorrected linear; + grad(p) Gauss skewCorrected linear; + grad(U) Gauss skewCorrected linear; +} + +divSchemes +{ + default Gauss skewCorrected linear; + div(phi,U) Gauss skewCorrected upwind; + div(phi,k) Gauss limitedLinear 1; + div(phi,epsilon) Gauss limitedLinear 1; + div(phi,R) Gauss limitedLinear 1; + div(R) Gauss linear; + div(phi,nuTilda) Gauss limitedLinear 1; + div((nuEff*dev(grad(U).T()))) Gauss linear; + div(U) Gauss skewCorrected lineaUpwind; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(nuEff,U) Gauss linear corrected; + laplacian((1|A(U)),p) Gauss linear corrected; + laplacian((voidfraction2|A(U)),p) Gauss linear corrected; + laplacian(DkEff,k) Gauss linear corrected; + laplacian(DepsilonEff,epsilon) Gauss linear corrected; + laplacian(DREff,R) Gauss linear corrected; + laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; + laplacian(phiIB) Gauss linear corrected; + laplacian(U) Gauss linear corrected; +} + +interpolationSchemes +{ + default skewCorrected linear; + interpolate(U) skewCorrected linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution new file mode 100644 index 00000000..b4c9e0ec --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/CFD/system/fvSolution @@ -0,0 +1,118 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 6 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + pFinal + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + U + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } + + k + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } + + epsilon + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } + + R + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } + + nuTilda + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } + + phiIB + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + "(U|k)Final" + { + $U; + tolerance 1e-05; + relTol 0; + } + + Us + { + $U; + tolerance 1e-05; + relTol 0; + } + +} + +PISO +{ + nCorrectors 2; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + +relaxationFactors +{ +/* + equations + { + ".*" 0.9; + } +*/ +} + +// ************************************************************************* // diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt new file mode 100644 index 00000000..55c1fc96 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/data/spheres.txt @@ -0,0 +1,10 @@ +1.58162822 1.14912017 0.00000000 1.95500000 +0.60412822 1.85931549 0.00000000 1.95500000 +-0.60412822 1.85931549 0.00000000 1.95500000 +-1.58162822 1.14912017 0.00000000 1.95500000 +-1.95500000 0.00000000 0.00000000 1.95500000 +-1.58162822 -1.14912017 0.00000000 1.95500000 +-0.60412822 -1.85931549 0.00000000 1.95500000 +0.60412822 -1.85931549 0.00000000 1.95500000 +1.58162822 -1.14912017 0.00000000 1.95500000 +1.95500000 -0.00000000 0.00000000 1.95500000 diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init new file mode 100644 index 00000000..c1fb1d62 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_init @@ -0,0 +1,113 @@ +# Insert a RBC into a box, then induce shear flow +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 300 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Sben equal ${Sn}*10 # analogous bending modulus of the parallel bond +variable Stor equal ${Sn}*10 # analogous torsional modulus of the parallel bond + + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary p f p + +# communication +newton off +communicate single vel yes + +# units +units micro + +#region and box creation +region domain block 0 45 0 22.5 0 22.5 units box +create_box 2 domain + +# neighbors +neighbor 1 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e03 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.5 +fix m3 all property/global coefficientRestitution peratomtypepair 2 1e-06 1e-03 1e-03 1e-06 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.1 0.1 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 2.65 0.0 1 1e50 1e50 + + +# particle template, insertion region and bond formation +fix pts1 all particletemplate/multiplespheres 15485863 atom_type 1 & +density constant 1098 nspheres 10 ntry 1000000 spheres file ../DEM/data/spheres.txt scale 1 & +bonded yes/explicit & +nbond_pairs 10 & +1 2 & +2 3 & +3 4 & +4 5 & +5 6 & +6 7 & +7 8 & +8 9 & +9 10 & +10 1 & +bond_type 1 + +fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0 + +region bc1 block 0 45 0 22.5 0 22.5 units box + +fix ins1 all insert/pack seed 32452867 distributiontemplate pdd1 & + maxattempt 10000 insert_every once overlapcheck yes vel constant 0.0 0.0 0.0 & + region bc1 particles_in_region 1 pos 22.5 11.25 11.25 ntry_mc 100000 + +run 1000 + +# set timestep +timestep 0.001 + +# check timestep +fix ts_check all check/timestep/gran 1000 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 10000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump*.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bondtype f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds*.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 100000 + +write_restart ../DEM/post/restart/liggghts.restart + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run new file mode 100644 index 00000000..822edd9e --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/in.liggghts_run @@ -0,0 +1,130 @@ +# RBC in a box sheared by fluid +log ../DEM/log.liggghts +thermo_log ../DEM/post/thermo.txt + +# bond parameter setup +variable rp equal 1.955 # sphere radius +variable ro equal 1 # bond outer radius multiplier +variable ri equal 0 # bond inner radius multiplier +variable rb equal 2*${ro}*${rp} # bond radius +variable lb equal 0.3090 # bond length multiplier +variable damp equal 25 # damping coefficient + +# material properties +variable G equal 0.005333 # in-plane shear modulus of the membrane +variable Y equal 0.020426 # in-plane Young's modulus of the membrane +variable Sn equal $Y/3.91 # analogous Youngs modulus of the parallel bond +variable St equal $G/3.91 # analogous shear modulus of the parallel bond +variable Ab equal (22/7)*${rb}${rb} # bond area +variable Sben equal ${Sn}/${Ab} # analogous bending modulus of the parallel bond +variable Stor equal ${St}/${Ab} # analogous torsional modulus of the parallel bond + +# atom type and style +atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6 +atom_modify map array + +# boundary +boundary p f p + +# communication +newton off +communicate single vel yes + +# units and processors +units micro +processors 1 1 2 + +# read the restart file +read_restart ../DEM/post/restart/liggghts.restart + +# neighbors +neighbor 1.0 bin +neigh_modify exclude molecule all +neigh_modify delay 0 + +# material properties +fix m1 all property/global youngsModulus peratomtype 7.5e02 7.5e04 +fix m2 all property/global poissonsRatio peratomtype 0.5 0.18 +fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.1 0.1 0.9 +fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.3 0.3 0.3 + +# pair style +pair_style gran model hertz tangential history +pair_coeff * * + +# bond parameters +bond_style gran +# N ro ri lb Sn_bond St_bond s_bend s_tor damp bn bt TYPE_OF_BOND sigma_break tau_break +bond_coeff 1 ${ro} ${ri} ${lb} ${Sn} ${St} ${Sben} ${Stor} ${damp} 3.0 1.0 1.0 1e50 1e50 + + +# reset timestep after reading the restart file +reset_timestep 0 +timestep 0.1 + +# cfd coupling +fix cfd all couple/cfd couple_every 10 mpi +fix cfd2 all couple/cfd/force + +# apply nve integration to all constituent spheres - single sphere treatment +fix integr all nve + +compute unwra all property/atom xu yu zu + +# variables for data processing +variable x1 equal c_unwra[1][1] +variable y1 equal c_unwra[1][2] +variable z1 equal c_unwra[1][3] +variable x2 equal c_unwra[2][1] +variable y2 equal c_unwra[2][2] +variable z2 equal c_unwra[2][3] +variable x3 equal c_unwra[3][1] +variable y3 equal c_unwra[3][2] +variable z3 equal c_unwra[3][3] +variable x4 equal c_unwra[4][1] +variable y4 equal c_unwra[4][2] +variable z4 equal c_unwra[4][3] +variable x5 equal c_unwra[5][1] +variable y5 equal c_unwra[5][2] +variable z5 equal c_unwra[5][3] +variable x6 equal c_unwra[6][1] +variable y6 equal c_unwra[6][2] +variable z6 equal c_unwra[6][3] +variable x7 equal c_unwra[7][1] +variable y7 equal c_unwra[7][2] +variable z7 equal c_unwra[7][3] +variable x8 equal c_unwra[8][1] +variable y8 equal c_unwra[8][2] +variable z8 equal c_unwra[8][3] +variable x9 equal c_unwra[9][1] +variable y9 equal c_unwra[9][2] +variable z9 equal c_unwra[9][3] +variable x10 equal c_unwra[10][1] +variable y10 equal c_unwra[10][2] +variable z10 equal c_unwra[10][3] +variable time equal step*dt + +fix extra1 all print 10 "${time} ${x1} ${y1} ${z1} ${x2} ${y2} ${z2} ${x3} ${y3} ${z3} ${x4} ${y4} ${z4} ${x5} ${y5} ${z5} ${x6} ${y6} ${z6} ${x7} ${y7} ${z7} ${x8} ${y8} ${z8} ${x9} ${y9} ${z9} ${x10} ${y10} ${z10}" & + file ../DEM/post/particle_position.txt screen no + +# check timestep +fix ts_check all check/timestep/gran 1 0.1 0.1 + +# screen output +thermo_style custom step atoms numbond +thermo 1000 +thermo_modify lost ignore norm no +compute_modify thermo_temp dynamic yes + +dump dmp all custom 1000 ../DEM/post/dump.liggghts & + id mol type x y z vx vy vz fx fy fz omegax omegay omegaz radius +compute b1 all property/local & + batom1x batom1y batom1z batom2x batom2y batom2z & + batom1 batom2 btype bforceX bforceY bforceZ +# x1 y1 z1 x2 y2 z2 id1 id2 bondtype f_bond[1] f_bond[2] f_bond[3] +dump bnd all local 1000 ../DEM/post/bonds.bond1 & + c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] & + c_b1[7] c_b1[8] c_b1[9] c_b1[10] c_b1[11] c_b1[12] + +run 1 + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/post/restart/.gitignore b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEM/post/restart/.gitignore new file mode 100644 index 00000000..e69de29b diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh new file mode 100755 index 00000000..94bd244f --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/DEMrun.sh @@ -0,0 +1,26 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# DEMrun script for redBloodCellShearFlow test case +# init redBloodCellShearFlow +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- include functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath="$casePath" +headerText="run_liggghts_init_DEM" +logfileName="log_$headerText" +solverName="in.liggghts_init" +debugMode="off" +#------------------------------------------------------------------------------ + +#- call function to run DEM case +DEMrun $logpath $logfileName $casePath $headerText $solverName $debugMode + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md new file mode 100644 index 00000000..94373848 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/README.md @@ -0,0 +1,9 @@ +## Single Red Blood Cell Dynamics in Shear Flow + +This case uses a continuous forcing immersed boundary solver on the CFD side and +a reduced order model of a red blood cell on the DEM side, cf. + +Balachandran Nair, A.N., Pirker, S. and Saeedipour, M., "Resolved CFD-DEM +simulation of blood flow with a reduced-order RBC model", Comp. Part. Mech. +(2021) + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh new file mode 100755 index 00000000..a0cd6f60 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/parCFDDEMrun.sh @@ -0,0 +1,59 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# parCFDDEMrun script for redBloodCellShearFlow test case +# run redBloodCellShearFlow CFD-DEM +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +. ~/.bashrc + +#- source CFDEM functions +source $CFDEM_PROJECT_DIR/etc/functions.sh + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" +logpath=$casePath +headerText="run_parallel_cfdemSolverIBContinuousForcing_redBloodCellShearFlow" +logfileName="log_$headerText" +solverName="cfdemSolverIBContinuousForcing" +nrProcs="2" +machineFileName="none" # yourMachinefileName | none +debugMode="off" # on | off| strict +testHarnessPath="$CFDEM_TEST_HARNESS_PATH" +runPython="false" +postproc="false" +#------------------------------------------------------------------------------ + +#- call function to run a parallel CFD-DEM case +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode + +if [ $runPython == "true" ] + then + cd $casePath + python results.py +fi + +if [ $postproc == "true" ] + then + #- get VTK data from liggghts dump file + cd $casePath/DEM/post + lpp dump* + + #- get VTK data from CFD sim + cd $casePath/CFD + reconstructParMesh -constant -mergeTol 1e-06 + reconstructPar -noLagrangian + foamToVTK + + #- start paraview + paraview + + #- keep terminal open (if started in new terminal) + echo "...press enter to clean up case" + echo "press Ctr+C to keep data" + read +fi + + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh new file mode 100755 index 00000000..a3efb45b --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/post.sh @@ -0,0 +1,31 @@ +#!/bin/bash +#------------------------------------------------------------------------------ +# post script for redBloodCellShearFlow test case +# run post-processing +# Achuth N. Balachandran Nair - October 2018 +#------------------------------------------------------------------------------ + +#- source CFDEM env vars +source ~/.bashrc + +#- source CFDEM functions +shopt -s expand_aliases +source $CFDEM_PROJECT_DIR/etc/bashrc + +#------------------------------------------------------------------------------ +#- define variables +casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")" + +#------------------------------------------------------------------------------ +#- get VTK data from liggghts dump file +cd $casePath/DEM/post +lpp dump* + +#- get VTK data from CFD sim +cd $casePath/CFD +reconstructParMesh -constant -mergeTol 1e-06 +reconstructPar -noLagrangian + +#- start paraview +paraview + diff --git a/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py new file mode 100755 index 00000000..f11fbfc1 --- /dev/null +++ b/tutorials/cfdemSolverIBContinuousForcing/redBloodCellShearFlow/results.py @@ -0,0 +1,66 @@ +#!/usr/bin/python + +import numpy as np +import matplotlib.pyplot as plt + +# Reading sphere data file +t,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,x9,y9,z9,x10,y10,z10 = np.genfromtxt('DEM/post/particle_position.txt', skip_header=0,usecols = (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30),unpack=True) +t = t*1e-06 + +# Calculating the centre of mass of RBC +xcm = (x1+x2+x3+x4+x5+x6+x7+x8+x9+x10)/10 +ycm = (y1+y2+y3+y4+y5+y6+y7+y8+y9+y10)/10 +zcm = (z1+z2+z3+z4+z5+z6+z7+z8+z9+z10)/10 + +# Calculating the gyration tensor +G11 = ((x1-xcm)*(x1-xcm) + (x2-xcm)*(x2-xcm) + (x3-xcm)*(x3-xcm) + (x4-xcm)*(x4-xcm) + (x5-xcm)*(x5-xcm) + (x6-xcm)*(x6-xcm) + (x7-xcm)*(x7-xcm) + (x8-xcm)*(x8-xcm) + (x9-xcm)*(x9-xcm) + (x10-xcm)*(x10-xcm))/10 +G12 = ((x1-xcm)*(y1-ycm) + (x2-xcm)*(y2-ycm) + (x3-xcm)*(y3-ycm) + (x4-xcm)*(y4-ycm) + (x5-xcm)*(y5-ycm) + (x6-xcm)*(y6-ycm) + (x7-xcm)*(y7-ycm) + (x8-xcm)*(y8-ycm) + (x9-xcm)*(y9-ycm) + (x10-xcm)*(y10-ycm))/10 +G13 = ((x1-xcm)*(z1-zcm) + (x2-xcm)*(z2-zcm) + (x3-xcm)*(z3-zcm) + (x4-xcm)*(z4-zcm) + (x5-xcm)*(z5-zcm) + (x6-xcm)*(z6-zcm) + (x7-xcm)*(z7-zcm) + (x8-xcm)*(z8-zcm) + (x9-xcm)*(z9-zcm) + (x10-xcm)*(z10-zcm))/10 +G21 = ((y1-ycm)*(x1-xcm) + (y2-ycm)*(x2-xcm) + (y3-ycm)*(x3-xcm) + (y4-ycm)*(x4-xcm) + (y5-ycm)*(x5-xcm) + (y6-ycm)*(x6-xcm) + (y7-ycm)*(x7-xcm) + (y8-ycm)*(x8-xcm) + (y9-ycm)*(x9-xcm) + (y10-ycm)*(x10-xcm))/10 +G22 = ((y1-ycm)*(y1-ycm) + (y2-ycm)*(y2-ycm) + (y3-ycm)*(y3-ycm) + (y4-ycm)*(y4-ycm) + (y5-ycm)*(y5-ycm) + (y6-ycm)*(y6-ycm) + (y7-ycm)*(y7-ycm) + (y8-ycm)*(y8-ycm) + (y9-ycm)*(y9-ycm) + (y10-ycm)*(y10-ycm))/10 +G23 = ((y1-ycm)*(z1-zcm) + (y2-ycm)*(z2-zcm) + (y3-ycm)*(z3-zcm) + (y4-ycm)*(z4-zcm) + (y5-ycm)*(z5-zcm) + (y6-ycm)*(z6-zcm) + (y7-ycm)*(z7-zcm) + (y8-ycm)*(z8-zcm) + (y9-ycm)*(z9-zcm) + (y10-ycm)*(z10-zcm))/10 +G31 = ((z1-zcm)*(x1-xcm) + (z2-zcm)*(x2-xcm) + (z3-zcm)*(x3-xcm) + (z4-zcm)*(x4-xcm) + (z5-zcm)*(x5-xcm) + (z6-zcm)*(x6-xcm) + (z7-zcm)*(x7-xcm) + (z8-zcm)*(x8-xcm) + (z9-zcm)*(x9-xcm) + (z10-zcm)*(x10-xcm))/10 +G32 = ((z1-zcm)*(y1-ycm) + (z2-zcm)*(y2-ycm) + (z3-zcm)*(y3-ycm) + (z4-zcm)*(y4-ycm) + (z5-zcm)*(y5-ycm) + (z6-zcm)*(y6-ycm) + (z7-zcm)*(y7-ycm) + (z8-zcm)*(y8-ycm) + (z9-zcm)*(y9-ycm) + (z10-zcm)*(y10-ycm))/10 +G33 = ((z1-zcm)*(z1-zcm) + (z2-zcm)*(z2-zcm) + (z3-zcm)*(z3-zcm) + (z4-zcm)*(z4-zcm) + (z5-zcm)*(z5-zcm) + (z6-zcm)*(z6-zcm) + (z7-zcm)*(z7-zcm) + (z8-zcm)*(z8-zcm) + (z9-zcm)*(z9-zcm) + (z10-zcm)*(z10-zcm))/10 + +# Calculating the eigen values of the gyration tensor +min_eig = [] + +for i in range(len(t)): + G = np.matrix([ [G11[i],G12[i],G13[i]], + [G21[i],G22[i],G23[i]], + [G31[i],G32[i],G33[i]] ]) + + eigenvalues = np.linalg.eigvals(G) + eig = np.min(eigenvalues) + min_eig.append(eig) + +plt.plot(t,min_eig,linewidth = 2) +plt.xlabel('Time (s)') +plt.ylabel('Min. eigen value of G') +plt.grid() +plt.show() + +plt.plot(z1,y1) +plt.plot(z2,y2) +plt.plot(z3,y3) +plt.plot(z4,y4) +plt.plot(z5,y5) +plt.plot(z6,y6) +plt.plot(z7,y7) +plt.plot(z8,y8) +plt.plot(z9,y9) +plt.plot(z10,y10) +plt.show() + +plt.plot(z1,x1) +plt.plot(z2,x2) +plt.plot(z3,x3) +plt.plot(z4,x4) +plt.plot(z5,x5) +plt.plot(z6,x6) +plt.plot(z7,x7) +plt.plot(z8,x8) +plt.plot(z9,x9) +plt.plot(z10,x10) +plt.show()