Tutorial: incompressible‐‐pimpleHFDIBFoam‐‐fallingParticleDistribution
This tutorial represents a situation where a number of non-spherical particles of different sizes is sedimenting in a rectangular domain.
Note that in order for the tutorial to be fast to evaluate even on personal computers, it is constructed as two-dimensional. However, the DEM part of the code is suitable only for three-dimensional simulations and the particles properties were adjusted in such a way that the tutorial gives plausible results.
The test domain is a hexahedron of dimensions (120mm x 0.1mm x 60mm) which has the Y-direction empty, that is, not active for the solution. The geometry is generated directly using the blockMesh OpenFOAM application and is displayed in the figure above, including its dimensions and boundary. As stated above, the front and back in the Y-directions are defined as type empty, the active boundaries treated as type wall are highlighted in green and preascribed with zeroGradient boundary condition for pressure and noSlip for fluid velocity. At the remaining type patch boundary, we fix the value of pressure, i.e., fixedValue set to uniform 0 is used, and prescribe a zeroGradient boundary condition for fluid velocity.
Details on the test geometry, mesh, and types of boundaries, see
"tutorialDirectory"/system/blockMeshDict
For details regarding boundary and initial conditions for the solved-for variables, see the files in the directory
"tutorialDirectory"/0.org/
The DEM solver is configured via the HFDIBDEMDict found at path
"tutorialDirectory"/constant/HFDIBDEMDict
The file used in this tutorial is:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: 8 |
| \ / A nd | Web: www.OpenFOAM.org |
| \/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object HFDIBDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
bodyNames ( "icoSphere" );
interpolationSchemes
{
U cellPointFace;
method line;
}
surfaceThreshold 1e-4;
stepDEM 0.02;
geometricD (1 -1 1);
recordSimulation true;
recordFirstTimeStep false;
nSolidsInDomain 1000;
outputSetup
{
basic true;
iB false;
DEM false;
addModel true;
parallelDEM false;
}
DEM
{
materials
{
particle
{
Y 1e7;
nu 0.3;
mu 0.75;
adhN 0;
eps 0.25;
}
wall
{
Y 1e7;
nu 0.2;
mu 0.75;
adhN 0;
eps 0.25;
}
}
LcCoeff 4.0;
collisionPatches
{
wall0
{
material wall;
nVec (-1.0 0.0 0.0);
planePoint (0.0 0.00 0.0);
}
wall1
{
material wall;
nVec (1.0 0.0 0.0);
planePoint (0.12 0.0 0.0);
}
wall2
{
material wall;
nVec (0.0 0.0 1.0);
planePoint (0.00 0.0 0.03);
}
wall3
{
material wall;
nVec (0.0 0.0 -1.0);
planePoint (0.0 0.0 -0.03);
}
}
}
virtualMesh
{
level 2;
charCellSize 0.0012;
}
icoSphere
{
fullyCoupledBody;
updateTorque true;
startSynced false;
rho rho [1 -3 0 0 0 0 0] 10000;
material particle;
U
{
BC noSlip;
}
bodyGeom convex;
interfaceSpan 1.0;
timesToSetStatic 80;
bodyAddition
{
addModel distribution;
distributionCoeffs
{
stlBaseSize 0.005;
addMode fieldBased;
fieldBasedCoeffs
{
fieldName lambda;
fieldValue 0.05;
}
addDomain boundBox;
boundBoxCoeffs
{
minBound (0 -0.001 -0.03);
maxBound (0.04 0.001 0.03);
}
scalingMode noScaling;
noScalingCoeffs{};
rotationMode noRotation;
noRotationCoeffs{};
}
}
}
// ************************************************************************* //
Going from top of the file, in the list bodyNames() the particle names that shall be active in the simulation are defined; if the particle is based on STL surface mesh, an STL file with a matching name has to be present in the directory constant/triSurface/. Here, we work with a particle named "icoSphere", which is linked to the file:
"tutorialDirectory"/constant/triSurface/icoSphere.stl
The global solver configurations are:
-
interpolationSchemessetting for the immersed boundary method -
surfaceThresholdthreshold for particle projection and interpolation schemes, -
stepDEMintegration step for DEM solver defined as fraction of the global CFD integration step. Consequently, the inverse value is the amount of DEM steps per one CFD iteration -
geometricDsolution directions, active directions: 1; inactive: -1 -
recordSimulationboolean option to separately record the position of particles at a given time -
recordFirstTimeStepboolean option to record initial position of particles added to the domain -
nSolidsInDomainmaximum number of particles that can be active within the domain, if not present, 1000 is used
The control of the outputs set up through the outputSetup dictionary where the following outputs may be enabled
-
basicsimulation time, and body velocities and locations per CFD step -
iBdetailed info regarding particle properties per DEM step -
DEMdetailed info regarding particle contact treatment -
addModeldetailed info regarding particle addition into the computational domain -
parallelDEMdetailed info regarding particle contact treatment from all subdomains for parallel computations
The DEM dictionary is used to set materials of solid phase and collision patches. The materials dictionary is split into sub-dictionaries where multiple materials might be defined using
-
Y- Young Modulus (material stiffness), -
nu- Poisson ratio, -
mu- static friction coefficient, -
adhN- normal adhesion coefficient, -
eps- coefficient of restitution (dissipation). Next the curvature coefficientLcCoeffrepresents the local curvature of the considered solids.collisionPatchesdictionary is split into sub-dictionaries while each sub-dictionary contacns a definition of a collision boundary for DEM which may or may not correspond to a system boundary. In this specific case, eachwallis defined as a dictionary consisting of -
materialenter the name of a material defined inmaterialsdictionary, -
nVecouter normal vector to the boundary, and -
planePointarbitrary point located in the collision boundary.
The virtualMesh dictionary is a setting of the contact treatment algorithm for the STL mesh-based solids. It is described by level - a decomposition level, similar to the corresponding snappyHexMesh setting, declaring how much the contact area will be refined; and by charCellSize, that is, the size of the characteristic computational cell for initial refinement of the contact area.
Finally, each particle listed in bodyNames() has to have its properties defined. This is done via a dictionary named according to the entry in bodyNames(), in this case, icoSphere{}. Within the dictionary, it is necessary to define: mode of particle motion, particle material and density, boundary condition to enforce on the fluid-solid interface, type of particle geometry, mode of particle addition into the domain, and additional settings for the solver numerics. Note that the combination of the particle listed bodyNames() and the corresponding dictionary acts as a template for generation and treatment of arbitrary number of particles based on the selected mode of addition.
In this tutorial, the particles motion is fully coupled with the fluid, and the particles have the density of 10000 kg/m3 and are of the particle material as defined in DEM.materials. The corresponding entries in the icoSphere template dictionary are:
-
fullyCoupledBody;if you wish to determine initial velocity you may enterfullyCoupledBody{velocit (0 1 0);}. - to enable particle rotation, define
updateTorqueand set ittrue - for particle to start synchronised with fluid velocity and rotation, set
startSynced truein this case we assume zero initial velocity for particle material particle-
rho rho [1 -3 0 0 0 0 0] *value*;, whererhois the standard OpenFOAMdimensionedScalarvariable.
The boundary condition at the fluid-solid interface is U{BC noSlip;}, which is the only value presently implemented.
From the point of geometry, the icoSphere particle is convex, which leads to the bodyGeom convex; entry. Also, there is an option to freeze the simulated particle at a place after a prolonged period of it not moving. This is done via timesToSetStatic 80, which means that the fullyCoupledBody will be converted to static after 80 DEM time steps of inactivity.
In this tutorial, a size distribution of particles of shape defined by /constant/triSurface/icoSphere.stl is generated at the during the run time of the simulation. The corresponding add model is addModel distribution and its settings are
bodyAddition
{
addModel distribution;
distributionCoeffs
{
stlBaseSize 0.005; //stating referential size of the STL file
addMode fieldBased; //selecting mode to condition particle addition
fieldBasedCoeffs
{
fieldName lambda; //name of the indicator field
fieldValue 0.05; //target integral of the field over the addDomain
}
addDomain boundBox; //choosing to create new particles within the bounding box
boundBoxCoeffs
{
minBound (0 -0.001 -0.03); //( 0 -1 -30) mm
maxBound (0.04 0.001 0.03); //(40 1 30) mm
}
scalingMode noScaling; //aded particles will not be additionally rescaled
noScalingCoeffs{};
rotationMode noRotation; //added particles will not be additionally rotated
noRotationCoeffs{};
}
}
This addModel in particular requires additional data in file
"tutorialDirectory"/constant/distributionDict
which contains:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: 8 |
| \ / A nd | Web: www.OpenFOAM.org |
| \/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object distributionDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.005; //conversion of particleSize to meters
distribution 7 //7 fractions, percentage representation
(
0
10
60
15
10
5
0
);
particleSize 7 //7 fractions, fraction sizes
(
0.4
0.8
1.2
1.6
2
3
4
);
The case setting is then concluded by defining the gravity in
"tutorialDirectory"/constant/g
with value (9.81 0 0 ); Also, as this tutorial is for CFD-DEM simulation the gravity for fluid flow has to be defined using fvOptions, see
"tutorialDirectory"/system/fvOptions
Case settings not explicitly mentioned above are common for all the OpenFOAM cases based on the pimpleFOAM solver.
The case is run using ./Allrun & or bash Allrun & script:
`#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/RunFunctions
rm -rf 0
cp -r 0.org 0
runApplication blockMesh # mesh generation, see system/blockMeshDict
application=`getApplication` # selects application (pimpleHFDIBFoam) from system/controlDict
runApplication $application # run the simulation itself
paraFOAM -touch # to create .OpenFOAM file for visualisation in paraview
after the simulation ends, you may proceed to visualization and optional post-processing.
Several approaches might be used to visualize the results of the CFD-DEM simulation depending on the required level of detail.
First is the minimalistic approach based only on Paraview and standard OpenFOAM-like postprocessing. Paraview with OpenFOAM plugin is started by running the command
paraFOAM
from the case directory, which results in opening of a Paraview with a loaded file, in this case: "fallingParticleDistribution.OpenFOAM." To set a suitable display, check the boxes marked within the red borders in the figure below and then click Apply.
Below, we show example of visualization of the simulation initial state. Note that U is the fluid velocity, and lambda is the indicator field marking the particles positions. The Transform fiter is used to display the velocity field and particles positions side by side. Also, for this test, it is suitable to rescale the velocity magnitudes as indicated in the Set Range box outlined in red.
Alternatively, the Threshold filter can be used to display both U and lambda fields in a single object as shown below. The Threshold filter is applied to the lambda field as shown in the box outlined in green.
- The more advanced approach enables the display of full particle geometries. However, it requires additional postprocessing steps before launching Paraview. First, the initial time step has to be removed by running
rm -rf 0/
in the case folder, to be followed by
python3 sync_time_levels.py
Next, to merge all individual particle STL meshes into one for each time level, the following script is to be executed,
python3 merge_STL_outputFiles.py
If everything proceeded correctly, the time levels present should be renamed to integer values starting with 0, and a new directory STLMerged/ should be present in the case folder.
With the data prepared, the results are displayed in Paraview similarly as above. The STL files for postprocessing STL_Results..stl have to be loaded from STLMerged/ as shown in yellow-outlined boxes below. The loaded STL files may be treated as common in Paraview while below, a possibility to change the solid particles color is depicted.
A usefull tool for vizualization and case debugging is highlighting the domain in which the particles are generated. This can be achieved using the Box filter as shown below.
Questions regarding the tutorial may be posed in Discussions here on github or via email openhfdib-dem@it.cas.cz.
openHFDIB-DEM wiki, pose questions in Discussions or via email.