Tutorial 4 - gnomeCreative/HYBIRD GitHub Wiki

logo

Granular Column Creation with DEM

This tutorial introduces the DEM part of the code. DEM is conceptually very different from LBM and is more difficult to initialize. The code needs to be given as input the initial position and velocity of every particle that it simulates. To do that, an additional file needs to be provided via the option particleFile. The structure of this file is described here. There are two ways to create this file: (i) using an external editor (a Python or MATLAB script, a spreadsheet...) or (ii) use a "recycle" file produced by another simulation. In this tutorial, we aim at the former strategy and produce an initial sample of particles using a MATLAB script.

setup_granular_collapse (a) creation of a dense granular column by compaction under gravity, and (b) collapse by removal of the right wall (Tutorial 5).

First, let's look at the geometry. We want to create a sample of densely packed particles of mean diameter $d = 0.002~\textrm{mm}$, starting from a loose collection of particles (a). The particles will have a radius randomly picked within an interval of $d \pm 5%$. The particles will be placed and released inside a container of lateral size $l = 0.035~\textrm{mm}$. The container width is slightly wider than the particles: $w = 1.2d$. This way, the particles will be arranged in a quasi-2D packing (b).

The script for particle generation is straightforward and can be found in the folder. It generates a file called particle_init.dat. The structure of this file is explained in detail in the pre-processing documentation. Once created, move this file into the main simulation folder.

We want to generate a dense sample, packed under the influence of gravity. To do this, we will use a simple linear contact model, whose stiffness does not refer to any specific material. The chosen value, $2 \cdot 10^{3} , \mathrm{N/m}$, is just high enough that the particles will not overlap by more than a few percent of their radius. We will exploit a numerical technique to obtain a dense sample: deactivating friction between particles. This will be done by setting to zero every parameter related to friction. In the next tutorial, friction will be reactivated to simulate a collapse problem.

Setting the Environment

All configuration files should be set inside the tutorials/granular_column folder. Here the DEM part is active, and we will need a particle file generated using the MATLAB script. There are no immobile particles, so the object file will be empty, as in the previous simulations (objectFile=./null.dat). Therefore, a null.dat file must still be provided. The configuration file is already ready for use (granular_column.cfg), but the most important features are described here. As an exercise, the user can try to reconstruct it from scratch.

We are using only the DEM solver, and we do not need the free-surface algorithm for this, so we set:

  • demSolver=1: DEM is on.
  • lbSolver=0: No LBM.
  • freeSurfaceSolver=0: This is related to the LBM part and is irrelevant here.
  • forceFieldSolver=1: We need an external forcing (vertical gravity).
  • staticFrictionSolver=0: We do not need static friction. In fact, turning this off will help to obtain a dense sample at the end of the simulation.

We need two types of output files: Paraview files for visualization, and recycle files that will be needed for restart (see Tutorial 5):

  • maxTime=1.0: The code will run for $1.0 \textrm{s}$ simulation time, which should be sufficient to obtain a dense sample.
  • screenExpTime=0.01: The code will show a state line on screen every $0.01 \textrm{s}$ simulation time.
  • partExpTime=0.01: The code will generate a particle Paraview file every $0.01 \textrm{s}$ simulation time.
  • partRecycleExpTime=0.1: The code will generate a particle recycle file every $0.1 \textrm{s}$ simulation time.

The boundary conditions are solid walls in every direction. This can be set by specifying:

  • boundary0 = stat_wall
  • boundary1 = stat_wall
  • boundary2 = stat_wall
  • boundary3 = stat_wall
  • boundary4 = stat_wall
  • boundary5 = stat_wall

The external forcing (gravity) is set by choosing:

  • forceX = 0.0
  • forceY = 0.0
  • forceZ = -9.806: This is gravity, acting in the $z$ coordinate.

Domain setup:

  • latticeSpacing = 0.00024: This is not an important parameter here, but it is good practice to choose a value much smaller than the minimum dimension of the simulation domain. In our case, it is ten times smaller than the dimension in $y$ (see below).
  • domainSizeX = 0.035: The domain in $x$ is $0.035 \mathrm{m}$.
  • domainSizeY = 0.0024: The domain in $y$ is $0.0024 \mathrm{m}$, which is $0.2$ times the mean particle diameter, so just a bit larger.
  • domainSizeZ = 0.2: The actual size of the domain in $z$ is irrelevant, as long as there is enough space for the particles. It is important to look at the initial particle file and check the maximum values of the coordinates in the $z$ direction to estimate the space needed.

The physical parameters are the following:

  • contactModel = LINEAR: We use a simple linear contact model (see the DEM documentation).
  • linearStiff = 0.2E-4: This is the contact stiffness in the normal direction $k^\text{n}_\text{L}$.
  • particleDensity = 2500.0: The particle mass density $\rho_\textrm{s}$.
  • restitution = 0.88: The restitution coefficient $\zeta$.
  • viscTang = 0: The damping ratio in the tangential direction $\alpha_\text{t}$.
  • frictionCoefPart = 0, frictionCoefWall = 0, frictionCoefObj = 0: We will deactivate friction to obtain the densest possible packing.
  • particleFile = ./particle_init.dat: The file containing the initial location and velocity of the particles, as well as their radius.
  • numVisc = 0.001: This is a stabilization parameter. Choose something just small enough so that it will not affect the results. To get an idea about the value, imagine the particles immersed in a fluid (e.g. air). This is the dynamic viscosity of that fluid.

The code is configured by setting the following options:

  • fluidTimeStep = 5.0e-5: Although LBM is not running, the code needs a value here, smaller than the output step screenExpTime, but larger than the DEM time step. Something of the order of $10^{-4}$ or $10^{-5}$ will usually work well.
  • criticalRatio = 0.002: This is the ratio between the DEM time step and the expected collision time. It should be kept small, but remember that the smaller it is, the longer the simulation. Something of the order of $0.01-0.05$ usually works well. This is the standard way in the code to set up the DEM time step. The alternative is to use the multiStep option, which is not covered in this tutorial.

Running the code and visualizing the results

To run the code, follow the procedure explained in the previous tutorials. The code will run and print some reference state variables on the screen. The goal here is to obtain a stable, dense packing, which is in static conditions. Therefore, we need to monitor the simulation to ensure that static conditions have been achieved. This can be done by looking at the information that the code prints during execution, see a sample in the figure below.

screen_granular_collapse A sample of the screen output of the DEM simulation of granular column creation.

There are two ways to check the state of the system. The first is to look at the maximum speed, named "MaxPSpeed". This is the velocity of the fastest particle. It is therefore a global parameter. If the simulation is run for long enough, this parameter should initially increase because the particles are collapsing from their loosest initial state to a dense packing. Once the particles have collapsed, they will collide dissipating energy, and therefore their velocity will gradually (but chaotically) decrease. Once "MaxPSpeed" is low, one can safely assume that the system is static.

A more rigorous way to look at the system state is to observe the value of the kinetic energy, which is also printed on screen as eKin. This is the overall kinetic energy of the particles. It does not suffer from the same chaotic behaviour observable with MaxPSpeed, and evolves more smoothly. Once it has stabilized to a very low value, the system is probably in a static state. However, there could still be one or a few particles moving, and therefore it is a good idea to also have a look at MaxPSpeed. The value of the kinetic energy is also printed on the file "energy.dat", which the code produces on the folder with the results. This file contains the time (first column) and many states parameters, including the kinetic energy (third column). One can plot the time evolution of the kinetic energy using this file, for example using Excel. The plot should look like this one:

energy_granular_column The evolution of the overall kinetic energy, over time.

The particles can be visualized as well, by loading the VTU Paraview files contained in the folder particleData. Once loaded, the spheres are automatically visualized as a collection of points, which are the centres of mass of the spheres. To give them the appearance of spheres, two options are available. The best results are obtained by loading a filter called Glyph, and by choosing as Glyph type the option Sphere, see blue marker. Many settings have to be adjusted for achieving a good result, most are highlighted in red. The options are Radius=1, Theta Resolution=10, Phi Resolution=10, Scalars=radius, Scale Mode=Scalar, Scale Factor=1, and Glyph Mode=All points. This will create well-resolved sphere with the right radius around every point. They can be used to visualize results, for example to illustrate the value of the particle velocity by colouring them using v or the other fields saved in the Paraview files.

sphere_visualization Particle velocities and their visualization using glyphs.

All particle files are very light, and can be easily used to create animations, either live by using the options in the Animation View panel of Paraview, or using the command File $\rightarrow$ Save Animation, which will create a video file (very useful for presentations or videos to be attached to online versions of papers).