# Predicting average void fraction and void fraction uncertainty in fixed beds of poly-lobed particles

Matthieu Rolland, Andriarimina Daniel Rakotonirina, Anaïs Devouassoux,

José Luis Barrios Goicetty, Jean-Yves Delenne, Anthony Wachs

► To cite this version:

Matthieu Rolland, Andriarimina Daniel Rakotonirina, Anaïs Devouassoux, José Luis Barrios Goicetty, Jean-Yves Delenne, et al.. Predicting average void fraction and void fraction uncertainty in fixed beds of poly-lobed particles. Industrial and engineering chemistry research, 2019, 58 (9), pp.3902-3911. 10.1021/acs.iecr.8b05557 . hal-02090443

HAL Id: hal-02090443

https://ifp.hal.science/hal-02090443

Submitted on 4 Apr 2019

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Predicting average void fraction and void
fraction uncertainty in fixed beds of
poly-lobed particles

Matthieu Rolland,\*,† Andriarimina Daniel Rakotonirina,‡ Anaïs Devouassoux,†
José Luis Barrios Goicetty,† Jean-Yves Delenne,¶ and Anthony Wachs§,||

†IFP Energies nouvelles, Process Intensification Department, Rond-point de l'Echangeur de Solaize, BP 3, 69360 Solaize, France.

‡IFP Energies nouvelles, Fluid Mechanics Department, Rond-point de l'Echangeur de Solaize, BP 3, 69360 Solaize, France.

¶IATE, INRA, CIRAD, Montpellier Supagro, Université Montpellier, 34000 Montpellier, France.

§Department of Chemical & Biological Engineering, University of British Columbia,
Vancouver, BC, Canada V6T 1Z3.

||Department of Mathematics, University of British Columbia, Vancouver, BC, Canada
V6T 1Z2.

*E-mail: matthieu.rolland@ifpen.fr

Abstract

Random packed beds of cylindrical, trilobic and quadrilobic particles in cylindrical
and bi-periodic containers are numerically studied using Grains 3D, a code based on
the Discrete Element Method (DEM) that resolves all inelastic collisions and simulates
dynamically the loading of packed beds. To mimic industrial or laboratory packing
procedures, particles initial position and orientation are random so that the same sim-
ulation repeated again yields a different packed bed structure and thus a different aver-
age void fraction. These "in silico" experiments aim at being able to optimize particle
shape in heterogeneous catalysis, in particular with respect to the corresponding bed
void fraction that is a critical parameter for pressure drop prediction. These "in silico"
experiments are deterministic and accurate but with differences due to the loading pro-
cedure. In this paper, we first present our assessment of the uncertainty on average void
fraction induced by (i) the initial random position and orientation of inserted particles
and (ii) the insertion zone size. Next we investigate the effect of particle shape, namely
cylindrical, trilobic and quadrilobic on the average void fraction as a function of particle
length and diameter, and of the container type. Simple correlations are proposed that
describe very well the simulations within the aforementioned uncertainty related to the
packing procedure. While the beds made with cylindrical particles are markedly denser,
the beds made of the trilobic and quadrilobic particles have statistically an identical
void fraction.

Introduction

Numerous chemical reactions are industrially performed using heterogeneous catalysis. Cat-
alysts pellets can be shaped as spheres or extruded shapes (extrudates) or molded shapes. ¹⁻⁵
Due to the use of extrusion machines, extrudates are cheaper to produce in high quan-
tities. They can have various shapes: cylinders, trilobes, and more recently quadrilobes.
Molded shapes often include holes to improve internal transport. The best catalyst shape
is a compromise between catalyst cost, catalyst efficiency, pressure drop, attrition, and bed
plugging. ¹⁻⁵ Thus, it is application dependent. The challenge to design a better shape is to
be able to predict the gains based only on the shape knowledge.

Catalyst effectiveness is a measure of internal mass transfer limitation. It is defined as
the actual reaction rate (in mol/s) divided by the reaction rate that would be achieved if
the concentration inside the pellet was homogeneous and equal to that of the surface. If the reaction is fast enough, reactants may be consumed faster than they diffuse in the pellet so that they have a lower concentration at the pellet center than at its boundary. The active (expensive) catalytic phase located at the pellet center is not used as efficiently as the one at the surface. The engineering pathways to improve effectiveness are: (i) improving effective diffusion in the pellet by changing the pore size distribution and (ii) changing the shape and/or size, and introducing holes, to reduce the volume to external surface ratio. For a given shape, the catalyst effectiveness can be numerically predicted by solving the diffusion equation in the grains assuming a known kinetic scheme like in Mariani et al.⁶ With a little less accuracy, it can be reasonably predicted for any particle shape without holes using the generalized Thiele modulus as proposed by Aris,⁷ that can be written for a 1st order reaction:

$$ \Phi = \frac{V_p}{S_p} \sqrt{\frac{k}{D_{eff}}} \qquad (1) $$

$$ \eta = \frac{I_1(2\Phi)}{\Phi I_0(2\Phi)} \qquad (2) $$

where $V_p, S_p, k$ and $D_{eff}$ denote particle volume, particle surface, intrinsic kinetic constant and effective diffusion coefficient, respectively. $I_n$ is the modified Bessel function of order $n$. Reducing the particle diameter results in an improvement of the catalyst efficiency due to a lower $V_p/S_p$, unfortunately at the cost of a higher pressure drop.

Gas-Liquid pressure drop prediction in reactors has been the subject of many publications. Pressure drop estimations are always performed using at some point the single phase predictions, so that for our purpose, optimizing trickle bed pressure drop is the same as optimizing single phase pressure drop (see for example Attou et al.⁸). Pressure drop predictions are usually performed using correlations with a form following the Ergun's equation:⁹
$$ \frac{\Delta P}{H} = \alpha \frac{\mu(1 - \varepsilon)^2 u}{\varepsilon^3 d_p^2} + \beta \frac{\rho(1 - \varepsilon) u^2}{\varepsilon^3 d_p} \quad (3) $$

In the formulation of equation 3, the pressure drop is the combination of a frictional viscous term proportional to the velocity and a term proportional to the velocity square accounting for flow direction and section changes.¹⁰ Many publications give values or formulas for the constants α and β depending on the particle shape, mostly based on experimental results⁹,¹¹,¹² and more recently numerical results.¹³ Nevertheless, there is so far no universal method to precisely predict the Ergun's equation coefficients based on particle shape only.

As it can be noticed in equation 3, the pressure drop presents a very strong dependency on the void fraction. The void fraction depends strongly on the loading procedure as well as on extrudates length distribution, leading to some discrepancies between operators and catalyst batches. As the differences between most efficient shapes are small, it is difficult to experimentally decouple shape and length effects when measuring the packed bed void fraction. Repetition effects are barely quantified and are usually neglected, although we have no information on their magnitude compared to differences between shapes. It would be very helpful to be able to predict the void fraction with confidence intervals for new shapes.

In hydrotreatment applications, industrial extruded catalyst pellets are typically 1 to 2 mm in diameter by 2-8 mm in length. Their shape is cylindrical, trilobic or quadrilobic with various lobe shapes. This rather small aspect ratio guarantee a good mechanical strength of the particles. The void fraction of packed beds depends on many particle features such as, e.g., elongation, angularity, slenderness, polydispersity and non-convexity. The effect of each parameter on the bed porosity is not easy to isolate experimentally or numerically. For this reason only simple shapes have been studied in detail in the literature. It has been shown for example that ellipsoids can be packed more densely than spheres.¹⁴ In the case of simple spherocylindrical rods (a cylinder with a hemisphere at each end) nematic transition (self
alignment of particles) is observed for different packing conditions.¹⁵ Theoretical models only exist for particles of aspect ratio close to one¹⁶,¹⁷ or for fibers corresponding to the asymptotic case of very large aspect ratio¹⁸ but no practical correlations are available for the lobed shapes and aspect ratio (2-4) of hydrotreatment extruded catalyst pellets.

Catalysts can be packed in two types of reactors: large industrial reactors and "tubes" reactors, differing only in their size, with significant differences though. Large industrial reactors are a few meters wide so that particles arrangement and void fraction hardly depend on wall-particle interactions. This case can be simulated in a bi-periodic geometry. Tubes reactors are much smaller with typical inner diameter less than 100 mm for multi-tubular industrial plants and even smaller for laboratory / pilot reactors. We will here focus on laboratory reactors whose diameters are most often in the range of 10 to 25 mm. Some reactors of diameter 2 mm are now commercialized by companies like Avantium (Flowrence technology). In this case, the local void fraction is strongly affected by the presence of walls and the void fraction radial profile and average value of course depend on the reactor inner diameter as well as on the particles shape. The presence of walls make the particles organize in circles with an orientation that tend to be perpendicular to the radius and more horizontal than vertical.¹⁹

To summarize, it is not yet possible to predict the void fraction (and the pressure drop) accurately enough to rank innovative catalyst shapes. Experiments do not allow to control all parameters so that drawing definite conclusions about the effect of catalyst shape is difficult. New and advanced numerical tools are required to optimize the particle shape “in silico” that allows a perfect control of all parameters. In this paper, we use DEM to estimate the void fraction for trilobic and quadrilobic shapes, as well as analyze our computed results to establish trends in void fraction dependence to particle shape, loading procedure and container type.
# Methods and material

## DEM with non-convex particles

Several numerical methods to produce a packing of spheres have been published. Thanks to its flexibility, the Discrete Element Method (DEM) can be extended to more complex shapes and thus will be utilized in this work. Other methods than DEM can be used (see for example$^{20,21}$). The DEM method$^{5,22-26}$ is a Lagrangian particle tracking method which computes the velocity, trajectory and orientation of each individual particle in the system. A key feature of any DEM tool is its ability to detect collisions, determine the contact point(s) and compute the resulting contact forces. Recent developments of DEM allow the use of non-spherical particles, such as the glued spheres model which is a loose approximation of a complex shape,$^{27}$ or by an accurate description of arbitrary convex particles.$^{24,26}$ Recent developments by our group$^{5,28}$ allows the simulation of non-convex particles modelled as a collection of convex particles. This method, called “glued convex”, is an extension of the glued sphere method of Nolan et al.$^{27}$ It allows the use of the existing methods, models and algorithms already implemented in our DEM code Grains3D$^{24}$ such as the equations of motion, contact forces, time integration, collision resolution and particularly the Gilbert-Johnson-Keerthi algorithm$^{29,30}$ for collision detection. Detailed information about our code Grains3D, its DEM and parallel computing features and the extension to non-convex shapes can be found in other publications by Wachs et al. and Rakotonirina et al.$^{5,24,25,28}$, including the contact force model. The relevant parameters are:

* the stiffness coefficient $k_n$ that controls the maximum amount of artificial overlap between 2 rigid bodies,

* the coefficient of restitution $e_n$,

* the Coulomb friction coefficient $\mu_c$.
Simulation principle

Fixed beds of non-convex particles are computed using Grains3D. An insertion volume (also called insertion zone or window) is defined at the top of the domain as illustrated in figure 1. The insertion window can be a box, a flat surface or a single point. The particles are inserted in the simulation in the following sequence:

* the code attempts to insert a single particle per discrete time,

* for the next particle to be inserted in the system, the code selects randomly its center's position in the insertion volume and its orientation,

* the next particle is actually inserted if it does not collide with any other already inserted particles; if it does, the insertion is unsuccessful,

* when an insertion is unsuccessful, the code tries to insert the same particle again over the next discrete times until the insertion is successful; at each discrete time, a new random position and orientation is selected,

* all particles are subjected to the gravity force and progressively leave the insertion volume.

A large insertion zone results in more particles inserted over a short period of time as the time step magnitude is very small in DEM simulations. During their free fall, the particles experience inelastic collisions with walls and other particles. The total kinetic energy of the system decreases exponentially with time once all particles have reached their pseudo-final position in the bed. Simulations are completed when the maximum of the particle velocity norm decays below a tolerance criterion chosen here as $10^{-5} m/s$. Simulation outputs of interest here are the final position and the final orientation of each particle. In the horizontal direction, the domain geometry can be either constrained with rigid walls or semi-infinite using bi-periodic conditions.
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of a box-like insertion window in DEM simulations.</figcaption></figure>

## Void fraction analysis

The average void fraction (porosity) is computed from the Grains3D output file by two methods: (i) performing a 3D discretization of the space and counting the number of cells occupied by particles. Provided sufficiently small grid cells, this method is very accurate but computationally expensive.²⁸ (ii) a simplified method that we will present now.

In a given control volume of diameter equal to that of the reactor, the volume occupied by the particles is equal to the number of particles in the volume times the volume of a single particle. Assuming all particles are fully included in the volume, void fraction is given by :

$$ \varepsilon = 1 - \frac{N V_p}{\Delta z S} \qquad (4) $$

where N, Δz, Vp and S denote total number of particles, height of the bed, particle volume and reactor cross-section respectively.
Our simplified method consists in substituting $\Delta z$ in equation 4 by the difference of the vertical positions of the center of mass of the highest and lowest particles in the volume. This approximation is valid when the number of the particles cropped at the extremities of the volume is small compared to the number of particles inside the volume, in other words when the control volume is long enough. Equation 4 then becomes :

$$ \varepsilon = 1 - \frac{V_p}{S} \frac{1}{s} \qquad (5) $$

where $s$ denote the slope of vertical position vs. particle rank line. So our method consists in sorting all the particles according to their center of mass vertical position $z$, plotting that vertical position against the particle ranking as illustrated in figure 2 and measuring the average slope. For a random packing, the ranking plot is almost a straight line. Non-linear trends in the ranking plot bring information about the structure: steps indicate “structured packing”; a changing slope indicates a change in the average void fraction. This method is as accurate as the discretization method when the control volume is large enough. In all simulations, the trend was linear indicating that the average void fraction is not sensitive to the falling distance that reduces as the bed is being constructed.

<figure><img src="image_3.png" /><figcaption>Figure 2: Example of plots resulting from the method (ii) of the void fraction analysis.</figcaption></figure>

A correct estimation of void fraction requires to discard a few layers (2-4) at the top and bottom of packing,¹⁹ avoiding end effects (flat bottom influence at the bottom and free surface at the top).
## Cases description

A first set of simulations is performed using bi-periodic boundary conditions as shown in figure 3(b). This simulates a semi-infinite container, and models the packing in a large reactor. The container size is set to 18 mm after checking that this parameter has no effect on the void fraction. Another set of simulation is ran in small size cylindrical reactors using solid walls as shown in figure 3(a). The vessel diameters are 14 mm, 16 mm and 19 mm.

<figure><img src="image_4.png" /><figcaption>Figure 3: Type of domains considered in this study: a) Packing of TL in cylindrical domain; b) Packing of TL in bi-periodic domain.</figcaption></figure>

Simulations are performed with the following shapes:

* cylinders with a circular cross-section, named cylinders and denoted CYL,

* cylinders with a trilobic cross-section, named trilobes and denoted TL. In our glued convex DEM, the trilobes are modelled by numerically gluing 3 cylinders and a prism of equilateral triangular cross-section,

* cylinders with a quadrilobic cross-section, named quadrilobes and denoted QL. In our glued convex DEM, the quadrilobes are modelled by numerically gluing 4 cylinders
and a prism of square cross-section.

The cylinders, trilobes and quadrilobes are illustrated in figure 4. The cross-sectional diameter of trilobes and quadrilobes is defined as that of the circumscribed circular cylinder (see figure 4(d)). For identical diameter and length, TL and QL occupy a volume of respectively 69% and 74% of the cylinder. The particles length is set to 3 mm, 4 mm or 5 mm and their diameter is varying in the range of 1.0 to 2.5 mm.

<figure><img src="image_5.png" /><figcaption>Figure 4: Particle shapes in this study.</figcaption></figure>

In each simulation at least 1000 particles are inserted. The parameters of all numerical
simulations are listed in table 1. The restitution coefficient $e_n$ has been computed from the
rebound height of spherical particles made of the same material as the extrudates (porous
alumina). The surface of the particle is very rough which explains the high value of the
friction coefficient $\mu_c$ set to 0.55. Please note that this value is not universal and depends
on the nature of the material of the two bodies that experience frictional contact. We set it
to 0.55 as a result of crude friction experiments conducted in our lab designed to determine
the onset of sliding of an extrudate set on a rough inclined metal plane. Critical angles for
onset of sliding were recorded between 26° and 32°, giving $\mu_c$ between $tan(26°)$ and $tan(32°)$
and an average value of about 0.55. This value is in line with measurements of 0.5-0.7 by
Sliney et al.³¹ Extrudates manufactured differently would give a different value. Based on
a conservative $v_0 = 2 \text{ m s}^{-1}$ maximum collision velocity, the stiffness coefficient $k_n$ is chosen
such that the theoretical maximum overlap of 2 rigid bodies is less than 0.5% of the radius of
the sphere of equivalent volume. The overlap value is chosen as a trade-off between volume
loses due to overlap and computing time. The value of 0.5% was chosen as the highest
value with an effect of void fraction less than 0.0001. With this value, the sharp angles
of the selected particles are very well accounted for. Similarly, the time step magnitude is
set about 20 times smaller than the theoretical collision duration (see Wachs et al.²⁴ for
more details about how to estimate DEM simulation parameters). An example of a packing
sequence is shown in figure 5 for 3 mm and 1.6 mm diameter TL packed in a 14 mm diameter
container.

Simulations with random insertion and data analysis

As mentioned earlier, the particles are inserted in the simulation with a random position and
orientation. Afterwards, the simulations and measurements are deterministic. Every packed
bed thus present a different void fraction due to different (random) initial conditions. As we
are interested in comparing the effects of shape on void fraction, we must be able to quantify
<table><caption>Table 1: Contact force model parameters and estimate of contact features at $v_0 = 2 \text{ m s}^{-1}$. $\delta_{max}$ denotes the theoretical maximum overlap of 2 rigid bodies and $R_e$ denotes the radius of a sphere of same volume as the particle.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$k_n$ (N m<sup>-1</sup>)</td><td>1 &times; 10<sup>5</sup></td></tr><tr><td>$e_n$</td><td>0.7</td></tr><tr><td>$\mu_c$</td><td>0.55</td></tr><tr><td>$\delta_{max}$ (m), $\delta_{max}/R_e$</td><td>1.5 &times; 10<sup>-5</sup>, 0.005</td></tr><tr><td>$T_C$ (s)</td><td>2.01 &times; 10<sup>-5</sup></td></tr><tr><td>$\Delta t$ (s)</td><td>1 &times; 10<sup>-6</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Packing sequence for 3 mm long x 1.6 mm diameter TL particles in a 14 mm diameter container at time t = 0.2, 0.5, 1 and 4.4 s. Colored by translation velocity norm.</figcaption></figure>
which part of the differences between two simulations are due to the shape or to the random insertion at the top of the domain.

## Repeating the packing

Several packings with the same set of 1000 particles are repeated for the 3 shapes (table 2). As the particle shape and dimension differ from one case to another, the average void fractions should not be compared for the moment but the reader should focus on the void fraction standard deviation ($\sigma < 0.0053$) which reads:

$$ \sigma = \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\varepsilon_i - \mu)^2}, \quad \text{where} \quad \mu = \frac{1}{N} \sum_{i=1}^{N} \varepsilon_i \qquad (6) $$

where $N$ and $\varepsilon_i$ stand for total number of simulations and void fraction of the simulation $i$ respectively.

Void fraction standard deviation accounts for both shape effects and packing repeatability. To estimate only the repeatability, we remove from all void fraction the average of each subset (shape effect) and compute the standard deviation of the whole ensemble (18 elements with an average of 0). We find the repeatability: $\sigma_1 = 0.0042$. At this point, it is worth reminding that once the particles are inserted in the simulation, the solver is deterministic: $\sigma_1$ is a measure of the effect of the random initial conditions.

<table><caption>Table 2: Repetition of random packing with identical particles.</caption><thead><tr><th rowspan="2">Shape</th><th rowspan="2">Length<br>[mm]</th><th rowspan="2">Diameter<br>[mm]</th><th rowspan="2">Number of<br>repetition</th><th colspan="4">Void fraction</th></tr><tr><th>Average</th><th>Std dev.</th><th>Min.</th><th>Max.</th></tr></thead><tbody><tr><td>CYL</td><td>3</td><td>1.6</td><td>7</td><td>0.3829</td><td>0.0053</td><td>0.3752</td><td>0.3921</td></tr><tr><td>TL</td><td>3</td><td>2.2</td><td>5</td><td>0.4127</td><td>0.0007</td><td>0.4118</td><td>0.4136</td></tr><tr><td>QL</td><td>3</td><td>2.2</td><td>6</td><td>0.4085</td><td>0.0050</td><td>0.4026</td><td>0.4150</td></tr></tbody></table>
Effect of insertion window size

We estimate the effects of the insertion window size for various geometric configurations of the container (cylindrical / bi-periodic and its size), and various particle shapes and sizes. In this work, we only use planar (2D) square insertion windows and an insertion point for reference (see figure 6). For 13 selected cases, we compare the void fraction computed using two insertion windows with different side length (see table 3). Simulations are not repeated.

<table><caption>Table 3: Effect of insertion window size on void fraction. Insertion window is 2D square of the length (0 mm to 10 mm). 0 indicates point insertion.</caption><thead><tr><th>Case</th><th>Window 1 size [mm]</th><th>Window 2 size [mm]</th><th>Void fraction difference (×10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>D14-CYL-L4, d<sub>p</sub> = 1.8</td><td>4</td><td>6</td><td>0.3</td></tr><tr><td>D14-CYL-L3, d<sub>p</sub> = 1.6</td><td>0</td><td>4</td><td>14.8</td></tr><tr><td>D14-CYL-L3, d<sub>p</sub> = 1.6</td><td>0</td><td>6</td><td>6.9</td></tr><tr><td>D14-CYL-L3, d<sub>p</sub> = 1.6</td><td>4</td><td>6</td><td>7.9</td></tr><tr><td>D14-CYL-L3, d<sub>p</sub> = 1.6</td><td>4</td><td>6</td><td>7.4</td></tr><tr><td>D16-CYL-L3, d<sub>p</sub> = 1.8</td><td>4</td><td>6</td><td>3.4</td></tr><tr><td>D16-CYL-L3, d<sub>p</sub> = 1.6</td><td>4</td><td>6</td><td>0.7</td></tr><tr><td>D16-CYL-L3, d<sub>p</sub> = 1.8</td><td>4</td><td>6</td><td>1.8</td></tr><tr><td>BIP16-CYL-L3, d<sub>p</sub> = 1.6</td><td>4</td><td>10</td><td>5.8</td></tr><tr><td>BIP18-CYL-L3, d<sub>p</sub> = 1.4</td><td>4</td><td>10</td><td>6.6</td></tr><tr><td>BIP16-CYL-L3, d<sub>p</sub> = 1</td><td>4</td><td>10</td><td>16.3</td></tr><tr><td>BIP18-QL-L3, d<sub>p</sub> = 1.9</td><td>7</td><td>10</td><td>1.9</td></tr><tr><td>BIP18-QL-L3, d<sub>p</sub> = 2.2</td><td>6</td><td>10</td><td>6.9</td></tr><tr><td></td><td></td><td>Average</td><td>6.2</td></tr><tr><td></td><td></td><td>Std dev.</td><td>4.9</td></tr></tbody></table>

According to an analysis of variance (ANOVA), the void fraction difference is statistically non zero. In average, a larger window results in a higher void fraction. We propose the following mechanism: a larger window results in more particles inserted simultaneously, leaving less time for a particle at the top of the stack to reach the most stable position before the arrival of the subsequent ones. The standard deviation on the void fraction difference is 0.0049. This value is rather small simply because the insertion window size variation is rather limited. In fact, this geometry has been chosen not to be too influential on the packing results.
<figure><img src="image_7.png" /><figcaption>Figure 6: Top view of two simulation domains with their corresponding insertion windows.</figcaption></figure>

Choosing the proper insertion geometry is a matter of compromise for several reasons. First, none of the methods is more realistic than another: in the laboratories, reactor loading is not standardized and is often manual. A change in particle size while keeping the same insertion window size results in a change in the number of particles that are inserted simultaneously, which yields more or less compact beds (the extreme being described as loose and dense packing in the literature). An obvious geometrical constraint is that the insertion window must be smaller than the reactor: smaller reactors need smaller insertion windows which leads to denser beds. This is similar to the reduction of the funnel diameter during an experimental loading. Last, a small insertion window requires a longer loading time, whereas a larger one permits a faster loading. In order to decrease the computing time and compute more cases, the simulations are performed with a medium size planar square insertion window (4mm and 6mm wide) that fits in all geometries. This choice will overestimate the void fraction compared to a point insertion and underestimate the void fraction for large particles. When looking at all the simulation data set, not only the cases in table 3, void fraction could not be correlated with, nor ranked against a ratio of box size to particle size. The insertion window size effect is a complex function of particle, reactor and box dimensions.
This effect is not random but appears so and we decided to model it as a Gaussian random variable. Its standard deviation $\sigma_2$ is equal to $1/\sqrt{2}$ of the standard deviation of the “void fraction difference” (see Appendix for details): $\sigma_2 = 0.00346 = 0.0049/\sqrt{2}$. $\sigma_2$ measures the unknown bias on the simulation induced by the choice of the insertion window size.

Simulated values for cylinders are on the dense side of the experimental values validating both the DEM simulation and our choices for insertion procedure. For an aspect ratio of 2, Zou and Yu's measurements³² can be extrapolated to 0.31 (dense) and 0.42 (loose) for large containers, to be compared with our simulation of 0.35 in the biperiodic domain. Similarly in packing of cylinders with an aspect ratio of 2 in tubes, they found void fraction range of 0.35-0.38 (dense) and 0.48-0.53 (loose) to compare with our simulations that give a void fraction near 0.39-0.40 for the same particle to tube diameter ratios.

## Overall uncertainty

An overall uncertainty on a single void fraction simulation result can now be estimated from $\sigma_1$ (random initial conditions) and $\sigma_2$ (bias induced by insertion window size). As both uncertainties are independent, a classical measurement statistic theory gives an estimate of the overall standard deviation: $\sigma = \sqrt{\sigma_1^2 + \sigma_2^2} = 0.0054$. An estimation of the overall uncertainty on a single measurement $I$ is $I = 2\sigma = 0.011$ (see Appendix for details). According to this analysis, there is a 95% probability that, given the output $\epsilon$ of a single experiment, the average void fraction of a large number of simulations falls in the interval $\epsilon \pm 0.011$ (with $\epsilon = 0.42$, this gives an estimate of 0.409 and 0.431). In other words, this corresponds to a relative uncertainty on the void fraction of less than 2.5%.
# Results

## Bi-periodic container

The average void fraction for various shapes, lengths and diameters simulated in a bi-periodic container is presented in figure 7. This case corresponds to large containers similar to industrial reactors. The void fraction is linearly correlated with particle aspect ratio $L_p/d_p$. Shorter, rounder particles pack more densely. Cylindrical particles present a lower void fraction and lower dependence on the aspect ratio than poly-lobed shapes. Surprisingly, the void fraction of trilobes and quadrilobes can not be distinguished.

<figure><img src="image_8.png" /><figcaption>Figure 7: Average void fraction in a bi-periodic container for particles of various shape, length and diameter. Dashed lines correspond to the incertitude computed earlier.</figcaption></figure>

Two correlations are proposed to predict the void fractions in large containers with an accuracy better than the overall uncertainty:

$$ \text{CYL: } \varepsilon = 0.289 + 0.033 \frac{L_p}{d_p} \qquad (7) $$

$$ \text{TL \& QL: } \varepsilon = 0.314 + 0.049 \frac{L_p}{d_p} \qquad (8) $$
In figure 7 the slope for the poly-lobed particles is much larger than that of the cylindrical particles. We suggest that during the packing, the lobes hinder rotation and result in a quick dampening of the vibrations induced by impacts. This results in less compact beds for poly-lobed particles.

Extending the trends to near spherical shape ($L_p/d_p = 1$) leads to a void fraction of 0.32 (CYL) and 0.36 (TL/QL) which are values close to dense packings of spheres. The ratio of the corresponding "asymptotic" solid volume fractions $1 - 0.32 = 0.68$ and $1 - 0.36 = 0.64$, respectively, is $0.64/0.68 \approx 0.94$, i.e. 94%. As expected, this ratio is higher than the ratio of volumes between cylinder and TL/QL equal to $\sim 70\%$, for the same diameter, as the lobes of one TL/QL particle can enter the "cylindrical envelop" of another TL/QL particle.

## Cylindrical container

### Cylindrical particles

The void fraction of a packed bed of cylindrical particles in a cylindrical reactor is in line with experimental measurements by Leva et al.³³ (our values are in the range of $d_p/D < 0.3$). The void fraction increases with the particle aspect ratio and seems to decrease with increasing reactor diameter $D$. However in the range studied, the effect of $D$ is barely larger than the repeatability. Following Leva et al.,³³ whose results suggest a proportional relationship to the inverse of vessel diameter, we propose the correlation in equation 9. It describes the whole data set with a maximum absolute error of 0.014 and a standard deviation of 0.006, which is about half of the uncertainty. Data points and equation 9 are plotted in figure 8(a). The correlation is as follows:

$$ \text{CYL: } \varepsilon = 0.315 + 0.0244 \frac{L_p}{d_p} + 0.141 \frac{L_p}{D} \qquad (9) \\ 14 < D[\text{mm}] < 19, \qquad 1 < L_p/d_p < 5, \qquad 3 < L_p[\text{mm}] < 4 $$
<figure><img src="image_9.png" /><figcaption>(a) Equation 9 vs. simulations. Dashed lines are parity ±I.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Comparison of experimental data with correlations</figcaption></figure>

Figure 8: Void fraction for packed beds of cylindrical particles in a cylindrical reactor for various reactor diameters, particle lengths and aspect ratios.

In our limited diameter range, a simplified correlation that does not take into account the reactor diameter D predicts the void fraction with only a slightly higher error and reads as follows:

$$ \text{CYL: } \varepsilon = 0.327 + 0.033 \frac{L_p}{d_p} \qquad (10) $$

$$ 14 < D[mm] < 19, \quad 1.67 < \frac{L_p}{d_p} < 2.86, \quad 3 < L_p[mm] < 4 $$

The maximum absolute error and standard deviation of equation 10 are 0.02 and 0.0077, respectively.

## Poly-lobed particles

The following linear correlation (equation 11) predicts the simulated void fraction with a lower accuracy than cylinders but equal to the overall uncertainty (see figure 9):
$$ \text{QL: } \varepsilon = 0.33 + 0.0328 \frac{L_p}{d_p} + 0.212 \frac{L_p}{D} \qquad (11) $$
$$ 14 < D[\text{mm}] < 19, \qquad 1.2 < \frac{L_p}{d_p} < 3.33, \qquad 3 < L_p[\text{mm}] < 4 $$

<figure><img src="image_11.png" /><figcaption>Figure 9: Void fraction for packed beds of quadralobal particles in a cylindrical reactor for various reactor diameters, particle lengths and aspect ratios: equation 11 vs. simulations. Dashed lines are parity $\pm I$.</figcaption></figure>

The results for TL particles are presented in figure 10. The following linear correlation (equation 12) also describes the data with an accuracy equal to the uncertainty:

$$ \text{TL: } \varepsilon = 0.345 + 0.0289 \frac{L_p}{d_p} + 0.15 \frac{L_p}{D} \qquad (12) $$
$$ 14 < D[\text{mm}] < 19, \qquad 1.2 < \frac{L_p}{d_p} < 3.3, \qquad 3 < L_p[\text{mm}] < 4 $$

A unified correlation predicting the void fraction for TL and QL regardless of the shape has the same accuracy as that of the TL (see figure 11). It is defined as follows (equation 13):
<figure><img src="image_12.png" /><figcaption>Figure 10: Void fraction for TL packed beds for various reactor diameters, particle lengths (not shown) and aspect ratio: equation 12 vs. simulations. Dashed lines are parity ±I.</figcaption></figure>

$$QL \ \& \ TL: \varepsilon = 0.329 + 0.0289 \frac{L_p}{d_p} + 0.15 \frac{L_p}{D} \qquad (13)$$

$$14 < D[mm] < 19, \quad 1.2 < \frac{L_p}{d_p} < 3.33,$$

$$2 < L_p[mm] < 4, \quad 1.2 < d_p[mm] < 2.48$$

# Discussion

## Effect of domain size in bi-periodic directions?

The bi-periodic simulation results presented so far have been obtained with a domain with a transverse size of 18 mm. 4 simulations have been repeated using smaller domains (8 mm and 10 mm) with CYL and QL with a particle aspect ratio of 3. The void fraction in smaller domains are within the repeatability of that in the large domain (transverse size of 18 mm). We have so far no indication of an effect of bi-periodic domain size in the range 8 mm to
<figure><img src="image_13.png" /><figcaption>Figure 11: Comparison of a unified correlation (equation 13) with numerical simulations for TL and QL in small cylindrical reactors.</figcaption></figure>

18 mm, although it would seem logical that packing structure and average void fraction will be affected at some point when decreasing the bi-periodic domain size. In presence of a stochastic injection process, it appears that performing simulations in the chosen domains does not impose any particular structure in the bed with a wavelength correlated to the transverse domain size. Simulation results would tend to indicate that even a transverse size of 8 mm is large enough to represent an infinitely large domain in the transverse direction. A more comprehensive study is required to conclude on this topic for example by studying the asymptotic convergence of void fraction when increasing bi-periodic domain size and compare it with the variability induced by the injection process.

## Remark on the effect of container size

For all three particle shapes (CYL, TL and QL), the void fraction is higher in small reactors than in semi-infinite vessels as expected. When the reactor diameter increases, none of the correlations presented in this work for cylindrical reactors converges to the correlation proposed for infinite vessels. This is not surprising as our cylindrical reactors are not large
compared to the particle length. In fact, the minimum $L_p/D$ in our simulations is 3/19 = 0.158, which is still high. To get asymptotically vanishing wall effects in a reactor, $L_p/D$ is probably required to be at least as small as 0.05. More simulations at larger reactor diameters and probably non-linear relationships would be necessary to propose a unified correlation.

## Effect of contact force parameters

Another important remark is that our data set is obtained for a single set of contact force pa- rameters, and in particular for single values of the restitution coefficient $e_n$ and the Coulomb friction force $\mu_c$. In beds of spheres, it is well documented that the void fraction increases with decreasing $e_n$ and increasing $\mu_c$. While the contact point and force chain network is undoubtedly different in a bed of cylinders, trilobic and quadrilobic particles compared to that in a bed of spheres, we anticipate the same qualitative trend. However, we have not performed any DEM simulations with varying $e_n$ and $\mu_c$ and cannot confirm this presumable qualitative trend for cylinders, trilobic and quadrilobic particles. At the very least, if this trend is indeed confirmed by DEM simulations, changing the value of $e_n$ and $\mu_c$ will change the value of the coefficients in the suggested correlations presented in equations 7, 8, 9, 10, 11, 12 and 13.

## Conclusion

DEM has been used to investigate packed beds of poly-lobed particles. Although our simu- lations are deterministic (and accurate), random input parameters (location and orientation of particles) as well as simulation parameters (insertion window) lead to an overall uncer- tainty on void fraction that has been estimated to be about ±0.011. A subsequent analysis of the void fraction and its dependence on the particle shape and reactor size showed that TL and QL present statistically identical void fractions. The effects of random insertion, i.e., filling procedure, in packed beds may mask the shape induced effect on void fraction
for optimized particles. We propose simple linear correlations to predict the void fraction for cylinders, trilobes and quadralobes in semi-infinite and small size cylindrical reactor that showed a good accuracy in the limited range of reactor diameter explored. More simulations and probably non-linear regressions are necessary to unify these correlations.

DEM simulations yield a rich set of information in the bulk of the bed. We have not yet exploited all these data. An extended analysis of the microstructure of the simulated bed would include extracting radial void fraction profiles in the case of cylindrical containers and orientation probability density functions for both bi-periodic domain and cylindrical containers, as in, e.g., Partopour et al.³⁴ and Dorai et al.¹⁹ The ability for 2 TL/QL static particles to exhibit multiple contact points due to their non-convex shape should also affect the force chain in the bed and may have an effect on the breakage of the catalyst particles. This is postponed to future work.

Ranking TL and QL and their chemical efficiency is not possible based only on void fraction. A precise knowledge of the relationship between shape and pressure drop is necessary to conclude. In this direction, another ongoing work in our group is to use particle resolved simulation (PRS) to compute the pressure drop in beds of poly-lobed particles. This is an extension of the work presented in Dorai et al.¹³ PRS has been shown to be easily extended to heat and mass transfer as the technique used to solve the momentum equations with solid obstacles, as, e.g., Immersed Boundary, Distributed Lagrange Multiplier/Fictitious Domain or lattice-Boltzmann, can be similarly used to solve other conservation equations.³⁴⁻³⁷ An integrated framework based on DEM and PRS for reactive flows, as demonstrated by Partopour et al.³⁴ and other groups, is currently at the final stage of development in our group⁵ and will enable us to assess the chemical efficiency uncertainty induced by the multiple input parameters: loading procedure, particle shape, contact force parameters, dimensionless momentum transfer numbers (Reynolds number), dimensionless mass transfer numbers (fluid/solid diffusivity ratio, Schmidt number, Damkohler number) and kinetic models.
# Supporting Information

## Appendix S1: standard deviation calculation

This information is available free of charge via the Internet at http://pubs.acs.org/.

# Acknowledgment

This work has been fully funded by IFPEN.

# References

(1) Moyse, B. Raschig ring HDS catalysts reduce pressure drop. *Oil Gas J.* **1984**, *82*.

(2) Cooper, B. H.; Donnis, B. B. L.; Moyse, B. Hydroprocessing conditions affect catalyst shape selection. *Oil Gas J.* **1986**, *84*.

(3) Afandizadeh, S.; Foumeny, E. Design of packed bed reactors: guides to catalyst shape, size, and loading selection. *Appl. Therm. Eng.* **2001**, *21*, 669-682.

(4) Mohammadzadeh, J. S.; Zamaniyan, A. Catalyst shape as a design parameter—Optimum shape for methane-steam reforming catalyst. *Chem. Eng. Res. Des.* **2002**, *80*, 383-391.

(5) Rakotonirina, A. Fluid-solid interaction in a non-convex granular media : application to rotating drums and packed bed reactors. Ph.D. thesis, University de Lyon, Lyon, France, 2016.

(6) Mariani, N. J.; Mocciaro, C.; Keegan, S. D.; Martínez, O. M.; Barreto, G. F. Evaluating the effectiveness factor from a 1D approximation fitted at high Thiele modulus: Spanning commercial pellet shapes with linear kinetics. *Chem. Eng. Sci.* **2009**, *64*, 2762-2766.
(7) Aris, R. On shape factors for irregular particles – I: The steady state problem. Diffusion and reaction. *Chem. Eng. Sci.* **1957**, *6*, 262–268.

(8) Attou, A.; Boyer, C.; Ferschneider, G. Modelling of the hydrodynamics of the cocurrent gas-liquid trickle flow through a trickle-bed reactor. *Chem. Eng. Sci.* **1999**, *54*, 785–802.

(9) Ergun, S. Fluid flow through packed columns. *Chem. Eng. Prog.* **1952**, *48*, 89–94.

(10) Larachi, F.; Hannaoui, R.; Horgue, P.; Augier, F.; Haroun, Y.; Youssef, S.; Rosenberg, E.; Prat, M.; Quintard, M. X-ray micro-tomography and pore network modeling of single-phase fixed-bed reactors. *Chem. Eng. J.* **2014**, *240*, 290–306.

(11) Carman, P. C. Fluid flow through granular beds. *Transactions-Institution of Chemical Engineers* **1937**, *15*, 150–166.

(12) Nemec, D.; Levec, J. Flow through packed bed reactors: 1. Single-phase flow. *Chem. Eng. Sci.* **2005**, *60*, 6947–6957.

(13) Dorai, F.; Teixeira, C. M.; Rolland, M.; Climent, E.; Marcoux, M.; Wachs, A. Fully resolved simulations of the flow through a packed bed of cylinders: Effect of size distribution. *Chem. Eng. Sci.* **2015**, *129*, 180–192.

(14) Donev, A.; Cisse, I.; Sachs, D.; Variano, E.; Stillinger, F.; Connelly, R.; Torquato, S.; Chaikin, P. M. Improving the density of jammed disordered packings using ellipsoids. *Science* **2004**, *990*.

(15) Ramaioli, M.; Pournin, L.; Liebling, T. M. Vertical ordering of rods under vertical vibration. *Phys. Rev. E* **2007**, *021304*.

(16) Man, W.; Donev, A.; Stillinger, F.; Sullivan, M.; Russel, W.; Heeger, D.; Inati, S.; Torquato, S.; Chaikin, P. M. Experiments on random packings of ellipsoids. *Phys. Rev. Lett.* **2005**, *198001*.
(17) Saint-Cyr, B.; Szarf, K.; Voivret, C.; Azéma, E.; Richefeu, V.; Delenne, J.-Y.; Combe, G.; Nouguier, C.; Villard, P.; Sornay, P.; Chaze, M.; Radjaï, F. Particle shape dependence in 2D granular media – CEGEO-shape Group. *Europhys. Lett.* **2012**, 98, 44008.

(18) Williams, S. R.; Philipse, A. P. Random packings of spheres and spherocylinders simulated by mechanical contraction. *Phys. Rev. E* **2003**, 78, 051301.

(19) Dorai, F.; Rolland, M.; Wachs, A.; Marcoux, M.; Climent, E. Packing fixed bed reactors with cylinders: influence of particle length distribution. *Procedia Eng.* **2012**, 42, 1335–1345.

(20) Lubachevsky, B. D.; Stillinger, F. H. Geometric properties of random disk packings. *J. Stat. Phys.* **1990**, 60, 561–583.

(21) Voivret, C.; Radjai, F.; Delenne, J.-Y.; El Youssoufi, M. S. Space-filling properties of polydisperse granular media. *Phys. Rev. E* **2007**, 76, 021301.

(22) Cundall, P. A.; Strack, O. D. A discrete numerical model for granular assemblies. *Geotechnique* **1979**, 29, 47–65.

(23) Cundall, P. A. Formulation of a three-dimensional distinct element model–Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts* **1988**, 25, 107–116.

(24) Wachs, A.; Girolami, L.; Vinay, G.; Ferrer, G. Grains3D, a flexible DEM approach for particles of arbitrary convex shape–Part I: Numerical model and validations. *Powder Technol.* **2012**, 224, 374–389.

(25) Rakotonirina, A.; Wachs, A. Grains3D, a flexible DEM approach for particles of arbi-
trary convex shape - Part II: Parallel implementation and scalable performance. *Powder Technol.* **2018**, *324*, 18–35.

(26) Seelen, L.; Padding, J.; Kuipers, J. A granular Discrete Element Method for arbitrary convex particle shapes: Method and packing generation. *Chem. Eng. Sci.* **2018**, *189*, 84–101.

(27) Nolan, G.; Kavanagh, P. Random packing of nonspherical particles. *Powder Technol.* **1995**, *84*, 199–205.

(28) Rakotonirina, A. D.; Delenne, J.-Y.; Radjai, F.; Wachs, A. Grains3D, a flexible DEM approach for particles of arbitrary convex shape - Part III: extension to non-convex particles modelled as glued convex particles. *Comput. Part. Mech.* **2018**, *5*, 1–30.

(29) Gilbert, E. G.; Johnson, D. W.; Keerthi, S. S. A fast procedure for computing the distance between complex objects in three-dimensional space. *IEEE J. Robot. Autom.* **1988**, *4*, 193–203.

(30) Gilbert, E. G.; Foo, C. Computing the distance between general convex objects in three-dimensional space. *IEEE Trans. Robot. Autom.* **1990**, *6*, 53–61.

(31) Sliney, H. E.; DellaCorte, C. Conference paper. NASA-TM-106348, DOE/NASA/50306-3, E-8121, NAS 1.15:106348 (Contract/Grant number: RTOP 505-63-5A). Prepared for U.S. Department of Energy, Conservation and Renewable Energy, Office of Vehicle and Engine R&D. October 1992.

(32) Zou, R.; Yu, A. Wall effect on the packing of cylindrical particles. *Chem. Eng. Sci.* **1996**, *51*, 1177–1180.

(33) Leva, M.; Grummer, M. Pressure drop through packed tubes. 3. Prediction of voids in packed tubes. *Chem. Eng. Prog.* **1947**, *43*, 713–718.
(34) Partopour, B.; Dixon, A. An integrated workflow for resolved-particle packed bed models with complex particle shapes. *Powder Technol.* **2017**, *322*, 258-272.

(35) Dorai, F.; Briquet, L.; Hammouti, A.; Rolland, M.; Wachs, A. Multi-Scale Simulation of Reactive Flow Through a Fixed Bed of Catalyst Particles. ASME 2014 4th Joint US-European Fluids Engineering Division Summer Meeting collocated with the ASME 2014 12th International Conference on Nanochannels, Microchannels, and Minichannels. 2014; p V01BT21A004.

(36) Deen, N.; Kuipers, J. Direct numerical simulation of fluid flow and mass transfer in dense fluid-particle systems. *Ind. Eng. Chem. Res.* **2013**, *52*, 11266-11274.

(37) Lu, J.; Das, S.; Peters, E.; Kuipers, J. Direct numerical simulation of fluid flow and mass transfer in dense fluid-particle systems with surface reactions. *Chem. Eng. Sci.* **2018**, *176*, 1-18.
List of Tables

<table><tr><td>1</td><td>Contact force model parameters and estimate of contact features at v<sub>0</sub> = 2 m s<sup>-1</sup>. &delta;<sub>max</sub> denotes the theoretical maximum overlap of 2 rigid bodies and R<sub>e</sub> denotes the radius of a sphere of same volume as the particle.</td><td>13</td></tr><tr><td>2</td><td>Repetition of random packing with identical particles.</td><td>14</td></tr><tr><td>3</td><td>Effect of insertion window size on void fraction. Insertion window is 2D square of the length (0 mm to 10 mm). 0 indicates point insertion.</td><td>15</td></tr></table>
List of Figures

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Illustration of a box-like insertion window in DEM simulations.</td>
      <td>8</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Example of plots resulting from the method (ii) of the void fraction analysis.</td>
      <td>9</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Type of domains considered in this study: a) Packing of TL in cylindrical domain; b) Packing of TL in bi-periodic domain.</td>
      <td>10</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Particle shapes in this study.</td>
      <td>11</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Packing sequence for 3 mm long x 1.6 mm diameter TL particles in a 14 mm diameter container at time t = 0.2, 0.5, 1 and 4.4 s. Colored by translation velocity norm.</td>
      <td>13</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Top view of two simulation domains with their corresponding insertion windows.</td>
      <td>16</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Average void fraction in a bi-periodic container for particles of various shape, length and diameter. Dashed lines correspond to the incertitude computed earlier.</td>
      <td>18</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Void fraction for packed beds of cylindrical particles in a cylindrical reactor for various reactor diameters, particle lengths and aspect ratios.</td>
      <td>20</td>
    </tr>
    <tr>
      <td>9</td>
      <td>Void fraction for packed beds of quadralobal particles in a cylindrical reactor for various reactor diameters, particle lengths and aspect ratios: equation 11 vs. simulations. Dashed lines are parity ±I.</td>
      <td>21</td>
    </tr>
    <tr>
      <td>10</td>
      <td>Void fraction for TL packed beds for various reactor diameters, particle lengths (not shown) and aspect ratio: equation 12 vs. simulations. Dashed lines are parity ±I.</td>
      <td>22</td>
    </tr>
    <tr>
      <td>11</td>
      <td>Comparison of a unified correlation (equation 13) with numerical simulations for TL and QL in small cylindrical reactors.</td>
      <td>23</td>
    </tr>
  </tbody>
</table>
For Table of Contents Only

<figure><img src="image_14.png" /></figure>