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

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 simulation repeated again yields a diﬀerent packed bed structure and thus a diﬀerent average 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 diﬀerences due to the loading procedure. In this paper, we ﬁrst 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 eﬀect 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. Catalysts pellets can be shaped as spheres or extruded shapes (extrudates) or molded shapes. [1][2][3][4][5] Due to the use of extrusion machines, extrudates are cheaper to produce in high quantities. 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. [1][2][3][4][5] 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 2 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. 6 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, 7 that can be written for a 1 st order reaction: where V p , S p , k and D ef f 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. 8 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. 10 Many publications give values or formulas for the constants α and β depending on the particle shape, mostly based on experimental results 9,11,12 and more recently numerical results. 13 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. 14 In the case of simple spherocylindrical rods (a cylinder with a hemisphere at each end) nematic transition (self 4 alignment of particles) is observed for different packing conditions. 15 Theoretical models only exist for particles of aspect ratio close to one 16,17 or for fibers corresponding to the asymptotic case of very large aspect ratio 18 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. 19 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.

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 • 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 µ c . 6

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.

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. 28 (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 : where N, ∆z, V p and S denote total number of particles, height of the bed, particle volume and reactor cross-section respectively. 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 : 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.

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. 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

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 δ 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. 13 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 (σ < 0.0053) which reads: where N and ε 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 repeatibility.
To estimate only the repeatibility, 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 repeatibility: σ 1 = 0.0042. At this point, it is worth reminding that once the particles are inserted in the simulation, the solver is deterministic: σ 1 is a measure of the effect of the random initial conditions.

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.  4.9 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.  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 (4 mm and 6 mm 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 σ 2 is equal to 1/

Overall uncertainty
An overall uncertainty on a single void fraction simulation result can now be estimated from

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.
TL & QL: ε = 0.314 + 0.049 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 polylobed particles.
Extending the trends to near spherical shape (L p /d p = 1) leads to a void fraction of 0.32 As expected, this ratio is higher than the ratio of volumes between cylinder and TL/QL equal to ∼ 70%, for the same diameter, as the lobes of one TL/QL particle can enter the "cylindrical envelop" of another TL/QL particle.

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. 33   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: 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): QL: ε = 0. 33   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: 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

Discussion
Effect of domain size in bi-periodic directions? 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 larged 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-

Conclusion
DEM has been used to investigate packed beds of poly-lobed particles. Although our simulations are deterministic (and accurate), random input parameters (location and orientation of particles) as well as simulation parameters (insertion window) lead to an overall uncertainty 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. 34 and Dorai et al. 19 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. 13 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. [34][35][36][37] An integrated framework based on DEM and PRS for reactive flows, as demonstrated by Partopour et al. 34 and other groups, is currently at the final stage of development in our group 5 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.