Rigid Body Dynamics Algorithm for Modeling Random Packing Structures of Nonspherical and Nonconvex Pellets.

Despite the common use of nonspherical catalyst pellets in chemical engineering applications, the packing structures of such pellets have not been as systematically studied and characterized as spherical packings. We propose a packing algorithm based on rigid body dynamics to simulate packing of nonspherical and possibly nonconvex pellets. The algorithm exerts a hard-body approach to model collision phenomena. The novelty is that the transition between moving and resting particles is controlled by a cutoff on the relative contact velocities, instead of artificially damping linear and angular velocities to stabilize the algorithm. The algorithm is used to synthesize packings of spheres, cylinders, and Raschig rings with tube-to-pellet diameter ratios 3-9.16. The packings are validated in terms of bulk porosity and radial void fraction distribution, finding satisfactory agreement with literature data. Denser packing structures are generated with high restitution coefficients and low friction coefficients. The confining tube walls play an important role, with highly fluctuating bulk porosities in narrow tubes.

Introduction

Fixed bed arrangements find extensive applications, particularly in reaction engineering, where they are employed as catalytic reactors for the transformation of reactants into desired products. The design of such systems is highly influenced by the structure of the packing matrix, which, in turn, is governed by the pellet and container size and shape, the loading method, and the subsequent treatment of the bed. Classical approaches in modeling fixed bed reactors, from the most simplistic model such as the pseudohomogenous plug flow model to the more improved Λ(r) model proposed by Winterberg and Tsotsas,[1] which imposes a Brinkman–Forcheimer-extended Darcy (BFD) model to account for the axial velocity field and to incorporate global and local bed properties, including bulk porosity and radial void fraction distribution. However, these pseudocontinuum models cannot provide an accurate prediction of the temperature field in tubular fixed beds, particularly for those with narrow to moderate tube-to-pellet diameter ratios (N), e.g., N ≤ 10, where the role of wall effects as well as pellet shape is completely neglected in such simplistic models.[2−4]

During the last decades, advanced numerical techniques such as computational fluid dynamics (CFD) and lattice Boltzmann (LB) methodologies are increasingly used to fully resolve the three-dimensional mass-, momentum-, and heat transport around (and inside) individual catalyst pellets inside a bed.[5−10] Before such detailed simulations can be executed, we need to know the positions and orientations of the individual catalyst pellets inside the bed. This has persuaded researchers to delve profoundly into topological details of random packing structures either through experimental techniques such as magnetic resonance imaging (MRI) and X-ray computer tomography (CT), e.g., refs (11−14), through numerical procedures in the form of in-house and ad hoc algorithms, e.g., refs (15−18), or even through commercial codes such as PFC3D, which is a commercial discrete element method (DEM) package, e.g., ref (59).

The majority of the prevailing efforts have concentrated on random packings of spheres, while application of catalyst pellets of nonspherical and often nonconvex shapes, such as cylinders, Raschig rings, pall rings, trilobes, etc., are becoming increasingly popular, particularly in chemical reaction applications because of their specific potential to enhance transport processes.[4] The amount of literature addressing structural properties of such nonspherical packings is scarce.[8,10,14,20,21] This can be ascribed to the cumbersome and complicated strategies necessary to predict the trajectories of nonspherical objects during the loading process, where the orientational freedom of such pellets may not only be very problematic in terms of collision modeling but also lead to exceptionally high computational expenses.[22,23]

The main aim of this contribution is to propose and examine a novel physics-based hard-body packing algorithm, capable of simulating the dynamics of the random packing process of nonspherical and even nonconvex pellets. The novelty of our method lies in the combination of an explicit treatment of instantaneous binary collisions with a global treatment of networks of contacts between multiple particles when they are reaching their resting state. The transition between moving and resting particles is controlled by a cutoff on the mutual contact velocity. The fidelity and robustness of the proposed packing procedure in reproducing the statistical mean properties of realistic fixed bed arrangements is then thoroughly investigated and validated using published experimental data. Before introducing the proposed algorithm, a detailed survey on previous research efforts in this field is given, highlighting their advantages and restrictions.

Literature Review

Topological features of random packing structures have been the subject of numerous analytical and experimental studies during the last five decades.[11,14,29−33,15,17,22,24−28] Nonetheless, only a minor selection of these studies has concentrated on evaluation and assessment of the spatial distribution of catalyst particles inside a bed, e.g., refs (11), (15), (31), and (33−35), whereas the amount of literature dealing with nonspherical pellets is scarce, e.g., refs (14), (20), and (33). Their generally opaque nature prohibits the exertion of conventional optical imaging techniques to address the spatial distribution of pellets of different shapes in random packing arrangements.[22] To circumvent this, several research groups have utilized advanced experimental techniques such as noninvasive imaging methods, e.g., MRI and 3d-CT.[12−14,29,36] For example, Sederman et al.[12] have used MRI in combination with image analysis techniques to characterize random packing structures of spheres with N = 9, 14, and 19. Ren et al.[13] have employed MRI, coupled with velocity encoding and pulsed-field-gradient nuclear magnetic resonance (PFG-NMR), to analyze the flow structure in random packings of cylinders and spheres with N in the range of 1.4 to 32. Zhang et al.[29] have coupled X-ray microtomography with a digital packing algorithm, which combines Monte Carlo and Distinct Element Methods, to reconstruct packing structures of equilateral cylinders with N = 12.8 from the microtomography images on a pellet-by-pellet basis. Baker et al.[14] have investigated the fidelity of an image-based meshing approach as a tool for reconstructing random packing structures through pellet-scale data extracted from noninvasive methods such as MRI and CT. Fundamentally, using such noninvasive methods, an automated image-based algorithm must be used to extract the pellet-scale information, which is challenging even for realistic packings of spheres (see refs (37) and (38)). However, for the case of nonspherical pellets, where their geometry embraces a variety of surface elements, i.e., flat and curved surfaces as well as corners, more complicated and rigorous reconstruction techniques are needed to accurately computate the positions and orientations of each pellet in a bed.[36]

On the whole, the complexities and costs associated with noninvasive experimental measurement of full 3d granular systems, together with the large computational costs for postprocessing, have persuaded researchers to seek alternative numerical methods. This has resulted in a multiplicity of packing algorithms, models and codes for generating random packing structures numerically.[11,15,17,31,34,35,39] Even though the specific details make each model unique, the majority of the prevailing algorithms can be classified into one of the following categories: (i) sequential deposition or deterministic algorithms, (ii) collective rearrangement algorithms, and (iii) physics-based methods.

Sequential deposition (SD) algorithms are principally initialized by either a sphere or sphere cluster. The packing structure is then generated on the basis of a procedure known as random settlement,[40,41] whereby the filling process is usually modeled by instantaneous placement of a new sphere in contact with either three spheres or with two spheres and the container wall, all of which are already fixed in their positions. A subcategory of deterministic approaches is the drop-and-roll type algorithm,[42] in which a packing is initially assembled in vertical direction and then is compacted through lateral displacements of pellets by imposing a gravitational force field. Coelho et al.[43] proposed a sequential packing model based on the successive deposition of grains in a gravitational field. Their model initiates with random placement of grains above the bed, and then the filling process is continued layer after a layer, until the grain particles reach a local minimum of their potential energy. The algorithm allows any displacement and rotation of particles that contributes to a lowering of their barycenters. Atmakidis and Kenig[44] have employed an improved version of this approach proposed by Kainourgiakis et al.[45] to generate random packing model of spheres with N = 5. Mueller[39] has investigated four different deterministic algorithms including modified Bennett,[46] layer, alternate, and percentage methods, to synthesize packing geometries of monosized spheres. What distinguishes these proposed approaches are the methods adopted to model the loading process, viz., the procedures devised for placement of new spheres in the alternative positions. Mueller demonstrated that the percentage model gives the best results for radial void fraction distribution out of the other approaches considered. The basis of this method is to locate a sphere at the lowest alternative vertical positions based on a prescribed percentage. It means that the resting site for a new sphere is determined by comparison between the percentages of accessible alternative sites and the prescribed one. Mueller suggested that the optimum value of this percentage corresponds to the lowest bed porosity; however, his further analyses demonstrated the severe dependence of this value on N, where for N ranging from 3.96 to 20.3, the optimum value of the percentage decreases from 70% to 10%. Furthermore, it was noticed that the accuracy of these deterministic approaches becomes lower as the bed-to-pellet diameter ratio increases. Mueller has improved his algorithm in another work,[34] by proposing a dimensionless packing parameter, thereby enhancing the procedure of the sphere’s sequential placement. Using this improvement, Mueller generated packings of monosized spheres with much better approximation of the radial void fraction profile, even for cases of large N-beds. Magnico[47] has used a modified version of the Bennet method, which was improved on the basis of the Mueller method[39] to account for wall effects, to generate random packings of spheres with N = 5.96 and 7.8. Basically, the distinct advantages of SD techniques are their low computational expenses as well as their intrinsic simplicities in terms of programming. However, as a major disadvantage, it is very difficult to manage the porosity of the resulting structures even with highly complex deposition algorithms.

Collective rearrangement (CR) algorithms form another category of packing simulation. CR is fundamentally regarded as a process whereby (i) the overlapping particles undergo a series of repetitive minor displacements so that all detected overlaps in the computational domain are removed (if the initial condition is overlapping) and (ii) the particles are migrated stochastically with the aim of decreasing the bulk porosity (if the initial condition is nonoverlapping). Liu and Thompson[48] have proposed a CR algorithm to generate random packing structures of spheres. The authors have investigated the influence of different boundary conditions on the generated packing structures. Maier et al.[49] have used a hard-sphere Monte Carlo algorithm to synthesize random packings of monosized spheres. In this algorithm, spheres are initially arranged at the intersections of a cubic lattice within a cylindrical tube with periodic boundaries set at its ends. Each sphere is then migrated stochastically, where a displacement is allowable if it does not lead to overlap with another sphere or the wall. The investigators benchmarked the radial void fraction profiles obtained from three generated-structures with N > 10 against experimental data. Sobolev and Amirjanove[50] have proposed a similar CR algorithm; however, their model benefits from a method allowing for a much denser initial arrangement of spheres in the container, thereby reducing the computational expenses significantly. The authors have shown that the density of simulated structures increases by repetition of the CR procedure, i.e., the number of packing trials. Freund et al.[5] have exercised a two-step Monte Carlo process based on the packing algorithm proposed by Soppe.[11] The CR algorithm exerted starts with initial placement of spheres inside a cylindrical tube, and then the packing is compressed through rearranging the spheres with an increased possibility into the direction of gravity. The packing process is stopped using convergence criteria that take into account changes in potential energy as well as mechanical stability.

A subcategory of CR algorithms is the digital packing approach, which is based on “pixelation” (2D) or “voxelation” (3D) of both objects and packing space.[20,31] The basis of a code called DigiPac, which has been developed and improved by Caulkin et al.,[20,31] is to allow random migration of particles, one grid for each time interval, on a cubic lattice. This random movement includes both directional and diffusive motions, imitating a random walk-based sedimentation model. In fact, the movement of particles over a grid facilitates the process of collision detection as it can be examined whether two particles occupy the same grid space at the same time. This interesting feature significantly reduces the computational expense of a typical run compared to other CR algorithms, whereby overlap detection is mostly undertaken by further mathematical analysis.[20] One of the most important advantages of this algorithm is to use digitization for representing packing objects, allowing for synthetic generation of random packings of even nonspherical particles. The investigators have inspected the validity of DigiPac by comparing the radial void fraction profiles, extracted from a number of simulated packing structures of both solid and hollow cylinders with 7 < N < 16, with their own experimental data. Baker and Tabor[51] have employed DigiPac to generate random packings of spheres including 160 particles with N = 7.14.

Fundamentally, CR algorithms have significant advantages compared to SD algorithms. For instance, the final bed porosity can be adjusted a priori, and a spatial correlation for particle size distribution as well as their initial placement can be prescribed. Nonetheless, the main disadvantage of such randomized particle packing algorithms is that they are neglecting the physical aspects of the packing process. Moreover, their computational demands are high due to the very slow process of convergence, particularly for very dense or very loose packings.[48]

The third classification, i.e., physics-based methods, embraces the whole procedure that realistically describes the interactions between pellets and between pellets and wall. The foundations of such procedures are rooted in Newtonian mechanics. Salvat et al.[15] have proposed a physics-based approach in the form of a soft-sphere algorithm to generate a packed bed structure of monosized spheres in a cylindrical container. The algorithm allows some interpenetration between particles but does not account for friction forces. The authors validated their model by comparing the predicted particle center distribution for N = 7.99 to experimental data by Mueller.[52] Dixon et al.[6] have employed an improved version of this algorithm using part of Mueller’s algorithm[39] concerning the initial placement of spheres at the base of the cylindrical container. The authors have then validated their packing generation procedure, for models of more than 1000 spheres, with N = 5.45 and 7.44, by comparing the predicted radial void fraction profile to the literate data reported by Benenati and Brosilow.[24] In a similar work, Behnam et al.[7] have validated their packing generation model using the experimental data reported by Mueller.[52] Siiria and Yliruusi[16] have developed a program based on Newton’s laws of motion to generate random packing of spheres. Their model accounts for all forces including gravity, collision forces, and friction.

The discrete element method (DEM) can be considered as a subcategory of physics-based methods. DEM (Cundall and Strack[53]) is conceptually related to molecular dynamics, in which the trajectories of individual particles are computed by evaluating all forces. Since the main emphasis of researchers in this area is to investigate dense granular flows, thereby generating the most packed particulate beds, most of investigators have thus implemented a time-driven (soft-sphere) class of DEM. Nowadays, DEM has become a common and reliable computational tool either to probe the dynamics of grains in a particulate bed[54] or to be coupled with CFD tools to investigate the hydrodynamics of a packed bed.[10,59,21,55,56] For example, Theuerkauf et al.[57] used DEM to investigate the local and global bed properties in narrow tubular fixed beds. The investigators have introduced DEM as a robust tool for generating random packing structures of spheres with moderate to low tube-to-pellet diameter ratios. Bai et al.[59] have exercised a commercial DEM code, PFC3D, to generate realistic random packings of spheres and cylinders with N ≤ 4 as a part of a CFD study of flow field and pressure drop in packed beds. The authors have highlighted the difficulties associated with simulation of nonspherical packing structures using DEM. Eppinger et al.[55] have exercised coupled DEM-CFD using a commercial CFD package, STAR-CD, to investigate the behavior of the flow field in DEM-generated packings of spheres with 3 ≤ N ≤ 10. The authors inspected the validity of their simulated models by comparing with published correlations concerning radial void fraction distribution. Yang et al.[58] have also employed the same DEM package as in the work of Bai et al.[59] to generate random packings of spheres with 3 < N ≤ 8, as a part of a CFD analysis of flow and thermal fields in packed beds. The investigators have shown satisfactory agreements between the radial void fraction data obtained from the models and published experimental data. In fact, DEM seems to be the most promising concept among the available methods for predicting both macroscopic and microscopic features of particulate flows. However, a particular challenge, not only in DEM but also in every numerical approach, is the accurate incorporation of nonsphericity, convexity, and nonconvexity of packing objects, as frequently encountered in practical situations. Only a few studies have dealt with nonspherical DEM, investigating various strategies to model nonspherical objects and contact detection algorithms. For a detailed review on the application of DEM in nonspherical particulate system, we refer to the works of Lu et al.[22] and Zhong et al.[23] These reviews have surveyed recent advances in 3d modeling of nonspherical objects, covering the main findings as well as the inherent difficulties concerned with nonspherical particulate systems using DEM.

The simplest nonspherical objects that have been incorporated into DEM are ellipses in 2D and ellipsoids in 3D, where a surrogate model of such particles can be simply represented by algebraic equations.[60,61] However, this method, i.e., modeling a nonspherical particle using algebraic equations, is restricted to particle shapes that can be expressed algebraically, as an (often relatively simple) equation. To date, one of the most frequent and straightforward approaches in nonspherical DEM is the so-called composite-sphere or glued-sphere method, in which the established framework of spherical DEM is applied to approximate nonspherical particles and their collisions during the packing process. A number of researchers have addressed the hindrances associated with the application of DEM glued-sphere approach such as its sensitivity to the parameters and the need for calibration.[62,63] Wu et al.[33] have employed this approach to investigate the behavior of cubical particle packings. The authors have investigated the influence of vibration conditions on packing densification. Dong et al.[10] have implemented the same approach using a commercial DEM code, STAR-CD (using the methodology proposed by Eppinger et al.[55]) to generate random packings of steatite rings (with do/d/h equal to 6.2/3.5/4.5 and 8/6/8 mm), as part of a CFD analysis of the flow field and heat transfer in pack beds. However, the main disadvantages of the composite-sphere method are (i) the computational expenses, which increase tremendously with the number of spheres required for 3d approximation of a nonspherical object (viz. the more realistic the geometry of a particle, the higher the number of spheres necessary per particle), and (ii) the occurrence of multiple contact points, which is inherent in nonspherical particles, which may lead to inaccurate force calculation.[22] This has urged several researchers to improve the glued-sphere DEM method.[63,64] Furthermore, alternative approaches to model collisional contact between particles have recently been proposed and applied to different problems, addressing different issues for modeling nonspherical particles.[22,65,66] Notably, the Gilbert–Johnson–Keerthi (GJK) distance algorithm is a fast computational procedure that can be used to detect overlap between nonspherical particles, which has recently been combined with a soft-particle approach for the contact model.[66] Of course, a GJK contact detection algorithm could also be combined with a rigid body dynamics approach (discussed next), but a major disadvantage of the GJK approach is that it is inherently limited to convex particles, i.e., particles without any holes and with only convex outer surfaces.

Despite the above seeds of research covering nonspherical pellets, there has been relatively little progress in this field. Several important modeling aspects, such as the contact detection procedure in such complicated systems, and the problem of force-torque coupling upon collision, require substantial research efforts. The main purpose of this paper is to address this complex problem, introducing a procedure of generating packings of nonspherical and even nonconvex particles. To this end, a physics-based packing algorithm based on rigid body dynamics (RBD) is developed and examined for generating random packings of catalyst pellets with different shapes, including spheres, equilateral solid cylinders, and Raschig rings. The essential features of our algorithm are (i) a realistic representation of nonspherical pellets using triangular face mesh, which allows for simulation of even nonconvex pellets with sharp edges, (ii) the exertion of a hard-body collisional model whereby avoiding the unphysically large overlap of particles that might be caused by artificially lowered spring stiffnesses, frequently employed in DEM simulations to prevent unfeasibly small time steps in the treatment of particle collisions, and (iii) an explicit way of modeling resting contacts between multiple particles based on relative velocities, which avoids artificial damping of linear and angular velocities. The fidelity of the approach in reproducing the topological properties of realistic packing arrangements is then discussed using comparisons between the predicted data obtained from computer-generated structures and published experimental data.

Numerical Setup and Model Formulations

Rigid Body Dynamics (RBD) and the Problem of Resting Contacts

RBD is frequently used in the field of graphic design and simulation.[67] Fundamentally, it is an analytical scheme capable of simulating the dynamic behavior of assemblies of arbitrarily shaped objects based on Newton’s laws of motion and Lagrangian mechanics. The background of RBD has been thoroughly investigated by several authors.[68,69] The main hypothesis of this approach is to consider the objects as rigid bodies, thereby facilitating the analysis of their motion by describing the translation and rotation of reference frames attached to each rigid body. Furthermore, contrary to the widely used time-driven DEM that benefits from the soft contact approach, RBD is more akin to hard contact methods, streamlining the process of collision analysis, particularly in a system of nonspherical objects, where any overlap between nonconvex objects makes the detection of contact points/edges very problematic.[22] RBD has been frequently used as an appropriate middleware framework in computer graphics and animation software such as FlipBook by DigiCel, Blender by the Blender Foundation, Maya by Autodesk, and Cinema 4D by Maxon. The application of RBD in the field of chemical reaction engineering has recently been introduced by Boccardo et al.,[8] where the authors used the open-source code Blender (which uses the Bullet Physics liberary) to synthesize packed beds of different pellet shapes such as spheres, cylinders, and trilobes. Following this, Partopour and Dixon[70] proposed an integrated workflow for resolved-particle fixed bed models with nonspherical pellet shapes. The authors have also used the Bullet Physics Library for generating packing structures of spheres, cylinders, Raschig rings, and quadrilobes with five holes. In the Blender software, a simplified approach based on damping the translational and rotational velocities of each object is used to speed up the process of getting objects to the resting position (see Blender reference manual at www.blender.org and Partopour and Dixon[70]). This simplistic approach is pursued to avoid overshoot and jiggling of objects in the resting contact condition which accordingly leads to a more stable simulation and convergence. However, the damping of linear and angular velocities is implemented on each “active” object during the simulation, i.e., at each time step, and causes a violation of the law of conservation of linear and angular momentum and energy over each time interval. In this paper we offer an alternative RBD-based packing algorithm in which no damping forces are applied to moving particles. Rather, the transition between moving and resting particles is controlled by a cutoff on the relative contact velocity followed by a detailed balance of constrained forces acting on each pellet to model the resting contact condition rigorously and to facilitate the stability of convergence in RBD simulations.

Description of a Nonspherical Pellet

In this work, a catalyst pellet is regarded as a nondeformable material, which is characterized by a translation vector x(t), indicative of the barycenter of a pellet in the world space, and a unit quaternion q(t) = [q0(t), q1(t), q2(t), q3(t)], with q02 + q12 + q22 + q32 = 1, which determines the orientation of a pellet around its center of mass in the world space. A rotation matrix R(t) that transforms the orientation of the particle from a body(-fixed) to the world space coordinate system can be expressed in terms of quaternions as

Using the quaternion-based approach, we can circumvent singularities (“gimbal lock”) inherent in the Euler angle method[71] and also diminish numerical drifts, and corresponding topological skewness, stemming from the direct exertion of rotational matrices in computing the orientation of a pellet during numerical simulation of the packing process. The other constant parameters describing the physio-mechanical properties of a catalyst pellet are its mass m, moment of inertia tensor Ib in the body space, collisional dissipation measured by a coefficient of restitution COR, and surface friction factor μ. Lastly, the most important part of the preprocessing step in RBD-based simulations is the 3d modeling of the pellet shape, accomplished by describing the particle boundaries in the body space. In order to model a surrogate for the pellets, a fast and simplified subdivision-based polygonal approach proposed by Loop[72] is adopted, by which the body surfaces of an object are approximated using triangular meshes. Having incorporated the general formula of packings such as sphere, equilateral solid cylinder, and Raschig ring into the Loop algorithm, the body space of surrogate geometries is modeled optimally and in the most smoothed fashion. Figure illustrates 3d models of the pellets investigated in this work.

Figure 1

3d models of catalyst pellets, represented by triangular face mesh: (a) sphere with dp = 10 mm, (b) equilateral solid cylinder with dp = h = 10 mm, and (c) Raschig ring with dpo/h/δ = 10/10/2 mm.

Rigid Body of Equation of Motion

The translational and rotational state of a moving object in the world space at time t is described here by a state vector, X(t), which is expressed aswhere x(t) is the position of the center of mass in world space, q(t) is the unit quaternion describing the object orientation, and P(t) and L(t) are linear and angular momentum, respectively. The equation of motion, EOM, for an assembly of rigid pellets would hence be the derivative of the state vector at time t:wherewhere i = 1, 2, ..., np denote the pellet index number, v(t) is the linear velocity of the center of mass of a pellet, ω(t) the angular velocity around its center of mass, F(t) is the total external force acting on pellet i, τ(t) is the torque on the object (around its center of mass), and I(t) is the moment of inertia tensor of pellet i, which can be simply computed from the body-fixed moment of inertia tensor Ib and the rotation matrix R(t). Note that both world space and body space (local) coordinate systems are required to rigorously describe the rotational motion of a nonspherical pellet in the world space during numerical packing simulations (see Figure a). Transformation between these two reference frames is conveniently performed using the rotation matrix at each time step. To compute the rotational motion of pellets in the world space system, the vector transformation approach is applied, which has been demonstrated to be more computationally efficient compared to the tensor transformation approach.[73]

Figure 2

Typical 3d-schematic of collision between two rigid cylinders: (a) spatial variables of colliding pellets and (b) normal and tangential forces acting on the contact point.

To track the change of the state variables of the catalyst pellets over time, the EOM for each pellet, which is conceptually an initial value ODE problem, is resolved using a midpoint scheme. The reason behind adopting such a method to handle EOMs is the order of precision of O(h3), which alleviates skewing effects in the orientations of pellets stemming from numerical drift during the packing process. The algorithmic update for the ith pellet can be expressed aswherewhere Δt is the time step at iteration m. To eliminate error accumulation in updating the orientation of a pellet, the quaternion of each pellet is renormalized after each time interval.[74]

Implementation of Force Fields and Torque

The total force acting on the ith pellet at time t, F(t) can be expressed aswhere W is a body force, in this case the weight acting on the center of mass of the pellet, causing it to accelerate at 9.81 ms–2 in the direction of gravity. This force is simply described aswhere ẑ is the unit vector in the direction of the z-axis in the world space coordinate system.

F(t) and F(t) describe the collision forces on pellet i due to its interaction with neighboring pellet or tube walls. In essence, for a typical oblique contact between two rigid objects, the contact force can be decomposed into normal and tangential directions. The normal contact force F(t) includes a repulsive force and viscous dissipation (associated with a certain coefficient of restitution). The tangential contact force F(t) is considered as a friction force. Figure shows a typical schematic of the forces acting on two colliding cylindrical pellets.

The total tangential force acting on the body of pellet i is expressed aswhere F(t) and F(t) are the tagential friction forces caused by lateral and bottom walls, respectively, and F(t) is the tangential friction force due to collisional contact between the ith and kth pellets. Fundamentally, the friction force between two objects resists the sliding motion of two contacted surfaces against each other. Here, the Coulomb friction model is employed to describe the friction force at contact time t. This model can be represented for colliding contact between pellets i and k asandwhere f(t) is the magnitude of the normal force at the contact point (discussed later), μd is the dynamic friction coefficient, t̂(t) is the unit vector in the world space coordinate system that friction force acts along, n̂(t) is the normal unit vector of the cell face or plane to which the contact point p belongs, and vp(t) and vp(t) are the surface velocities at the contact points on the ith and kth pellets, respectively, which can be mathematically expressed by

In order to evaluate if dynamic friction should occur at the contact point, a threshold εf, based on the relative tangential velocity of the colliding pellets on the face cell or plane, i.e., |vp·t̂|, is considered, thereby restricting the role of this force in the analysis when packings are being stabilized. It is worth remarking that both F(t) and F(t) can be computed using a simpler procedure by setting vpL = vpB = 0. Using the scalar vp(t) = vp(t)·n̂(t), which describes the relative contact velocity in the direction of n̂, we can classify the behavior of collision, into three categories. Fundamentally, positive values of vp(t), i.e., vp(t) > εp, mean that the pellets are separating and the contact point is vanishing after tc. In this condition, the ODE solver can still continue the computational procedure. A zero value within a numerical threshold, −εp < vp(t) < εp, indicates that the pellets are neither approaching nor receding, which describes a situation called resting contact. This is a complicated condition for modeling support and counter forces at contact points (this is meticulously tackled in section ). The last possibility, i.e., a negative value of this quantity, vp(t) < –εp, corresponds to the case of colliding contact. These alternative cases are illustrated in Figure .

Figure 3

Three alternatives for collision phenomena in a 2d vertex/face contact:(a) pellets are separating, (b) pellets are in resting contact, and (c) pellets are in colliding contact situation.

Similarly, the total normal force acting on the body of ith pellet, i.e., F(t), is described by the following relation:where the first two terms in the right-hand side of eq describe the normal force acting on the contact points of the ith pellet while in collisional contact with lateral and bottom walls, respectively, and the rightmost term, i.e., F(t), is the normal force due to interaction between the ith and kth pellets. Conceptually, F(t) and F(t) are always directed at the normal of the confining walls at the contact points and only affect the pellets, while F(t) acts on the contact points in the direction of normal vector of either the cell face or the plane passing through the two crossed cell edges, whereas two pellets are colliding (the procedure of contact time/point evaluation is explained in section ).

Conceptually, upon a collision, the normal relative velocities of the colliding pellets are suddenly reversed in a discontinuous way so that the bodies do not interpenetrate. Such a situation would violate the assumption under which the ODE solver is able to work, namely, a smooth variation of X(t) over the time intervals. Therefore, the ODE solver for the pellets engaged in collision conditions, whether it is colliding or resting contact, ought to be stopped, and another computational procedure should be pursued to trace the trajectories of the pellets. We note that this computational procedure is founded on the state vector of colliding objects at a time step before the contact time tc, implicating that we assume that (i) the state vectors of colliding objects do not change during the collision period, i.e., tc ϵ [t, t], and that (ii) the state vectors of objects in collisional contact change discontinuously at the end of the collision period, i.e., at t.

In this study, we will work with the impulse due to a colliding contact when vp(t) < – εp. For such a collision, taken place at time tc within the interval [t, t], with the net impulse is defined as

The above equation is nothing but the statement that the total impulse of all collisional forces on pellet i corresponds to the difference in linear momentum of pellet i. As for the case of friction force, consider that two pellets i and k are in colliding contact at the vertex point p in their body spaces at time tc, where p(tc) = p(tc) = p. In this condition, the velocities of the collided objects need to instantly undergo a drastic change to prevent the bodies from interpenetration. The impulse, J(tc), for this case can be expressed by the following relation:andwhere j and j are undetermined scalars representing the magnitude of impulse in normal and tangential directions with respect to the contact face, respectively. It is worth remarking that if the impulse of J acts on pellet i, the body of pellet k would then be subjected to an equal but opposite impulse of −J to satisfy the momentum conservation law. The impulse acting on pellet i (due to its collision with pellet k) produces an impulsive torque, τ, which is expressed as

In order to evaluate the magnitude of the impulse in the normal direction, j, an empirical model describing the colliding contact is used, where the relative normal velocities of the contact vertices before and after collision are connected by the coefficient of restitution, COR. This empirical law can be expressed aswhere vp– and vp+ are the relative velocities at the location of the contact in the n^(tc) direction before and after the collision time tc ϵ [t, t], respectively. The coefficient of restitution lies between 0 and 1. When it is equal to 1, the collision is perfectly elastic and no kinetic energy is lost, while COR = 0 results in vp+ = 0, corresponding to a plastic impact, where maximum kinetic energy is lost. Substitution of vp– and vp+ in eq , together with a convenient mathematical procedure (see Appendix A in the Supporting Information for a derivation), results in the following equation for the magnitude of normal impulse j:andwhere I0 is the unit matrix and K is a collision matrix for the colliding contact between pellets i and k. In case pellet i is in colliding contact with a confining wall, whether lateral or bottom wall, the magnitude of impulse can be derived in a straightforward way. For example, j, which describes F(tc), can be expressed bywhere

Suppose that the relative tangential velocity at the contact point p, where the ith and kth pellets are in colliding contact, vp–(t) = vp–(t). t^(t) is nonzero. The magnitude of the tangential impulse j acting at contact time tc ϵ [t, t] is defined by the product of the dynamic friction coefficient μd and j using Coulomb’s friction model as

However, if vp is very small, j calculated by eq may reverse the sign of vp+(t), resulting in unphysical motion. To circumvent such an incidence, we assume that j = 0 if vp–vp+ ≤ 0.

Having computed the total impulse, J(tc), for a typical colliding contact between the ith and kth pellets, the discontinuous changes of linear and angular velocities of the pellets after collision can be expressed asThe same relations can also be written for the linear and angular velocity differences of pellet k by substituting −J(tc) in the above equations.

Problem of Resting Contact

One of the most problematic conditions occurring during the numerical simulation of the packing process is the resting contact, in which bodies are neither colliding nor detaching at their contact points, i.e., −εp < vp(t) < εp. This situation may cause a divergent number of collisions in a finite time span, which hence needs a specific procedure to handle the contacts, such that the pellets are kept fixed in their positions without any interpenetration. Consider a configuration of n contact points at which bodies are in resting contact (see Figure ). We assume that in rest there are no tangential forces, viz., static friction is omitted, and therefore the resting contacts can be resolved by specifying a set of forces acting normal to the contact surface, i.e., φn^(tc), chosen in such a way that they obey three constraints: (i) they should prevent the pellets from interpenetration, (ii) they should act repulsively in a way that holds the bodies together in contact, and (iii) they should become zero at the moment the bodies begin to detach again.

Figure 4

Typical schematic of resting contact condition.

It is worth noting that all φ need to be determined simultaneously, since a force acting on one contact point may influence objects involved on other contact points. To this end, we introduce a distance function, dp,(t), describing the distance between the resting contact points from two pellets i and k in the collision period. Considering the notations exercised before in describing colliding contacts, the distance function for a typical resting contact point p (whether it is a vertex/face or edge/edge contact), at the contact time tc ϵ [t, t], can be expressed by dp,(tc):

The first constraint is to prevent bodies from interpenetration. Since both dp,(tc) and ḋp,(tc), where the latter describes the relative contact velocity in the direction of n^, vp(tc), are intrinsically zero (within the numerical threshold), we need to ensure that the second derivative of the distance function, measuring how two pellet surfaces accelerate toward each other at the contact point, to be equal to or greater than zero (within a numerical tolerance), i.e., dp,(tc) ≥ 0. The expression, describing dp,(tc) can be simply derived by two times differentiating eq :andwhere v̇p(tc) and v̇p(tc) are the accelerations of the contact points in the world space coordinate system.

It is worth remarking that if a pellet is in resting contact with the tube wall, eq can be simplified by setting the rightmost term to zero, viz., n̂(tc) = n(tc) = 0.

On the basis of the second constraint, each of the contact forces needs to act outwardly which means that φ(tc) ≥ 0. The third constraint can simply be described in terms of φ and dp,(tc). Since the contact force must possess a value of zero if the resting contact starts breaking, it implies that φ(tc) must be zero when a resting contact is broken. This constraint can be expressed by the following relation:The magnitudes of φ(tc) can thus be computed such that the mentioned constraints are satisfied for each contact point. Suppose that for a configuration of n (resting) contact points at tc ϵ [t, t], the distance acceleration for each of the contact points depends on the normal forces acting on all contact points together with some other constant forces and velocity-dependent terms related to the objects involved in that particular contact point at the precollision time step, t. In that case, we can express dp,(tc) in terms of all unknown constraint forces, including those directly affecting the pellets i and k aswhere a describes how a unit change in φ(tc), if pellet i is in direct contact with pellet j as well, can affect the contact point p. A more detailed derivation is given in Appendix B of the Supporting Information.

If we define Dp as the column vector of dp,(tc) for n resting contact points, at time tc, and as the column vector of all constraint forces φ(tc), then the problem of resting contact for n contact cases can be described in matrix form aswith the following constraints:with the understanding that eqs and 34 apply to each component of the column vectors. Equations –35 can be regarded as a quadratic programing problem, QP, in which we attempt to minimize the quadratic term ψTDp subject to the conditions ψ ≥ 0 and Dp ≥ 0. In this work, the procedures proposed by Baraff[75] is used to solve the problem, which has been proven as a reliable and fast solution algorithm.

Collision/Contact Detection Philosophy

The collision/contact detection procedure exerted here is based upon scrutinizing the proximities between body spaces of approaching rigid objects. This is implemented through a coherence-based two-layer search procedure, whereby the proximity between all pairs of face meshes are assessed to see if the minimum distance between their bounding volumes are lower than a preset threshold. The first layer of the search algorithm employs a series of criteria to narrow down the list of pairwise collision/contact possibilities at each time interval. These criteria are as follows:

For spheres:

For cylinders and Rasching rings (L = 2R):where λ = 0.1Rp is a proximity threshold considered between two approaching objects and Δ|| is the change of distance between barycenters of two objects in two successive time steps. The result of this search layer, in the form of a sorted list of pairs of pellets that are possibly colliding, is then fed to a supplementary search engine, thereby determining the time of contact as well as details of intersections between the face meshes of collided objects. Henceforward, to avoid any probable interpenetrations of the pairs of pellets reported by the first search layer, the midpoint scheme is run using a much smaller time interval, δt. To estimate δt, let us consider two pellets i and k, detected by the first collision detection search procedure, with a maximum relative surface velocity of |vmax(t)|. For the maximum allowable proximity of 0.2λ, the second collision search subalgorithm can initiate at

Using this time interval hereafter, the inner search layer is executed in the form of a sweep/sort algorithm, by which a sorted list of data including the points and cell face with minimum distance to the body spaces of each pair of listed objects is created. The list is then scanned for either an overlap or a contact between two meshes, viz. vertex/cell face or edge/edge contacts. A simple penalty method based on the distance between the reported contact points and the center of mass of the pair of pellets in contact is then exercised to avoid any interpenetration. Following this, the algorithm returns to the latest precollision time, allowing the postcollision calculation to proceed by computing the collision responses based on the type of collision. Figure illustrates how the two-layer collision detection procedure works.

Figure 5

Typical representation for two-layer collision search procedure.

Packing Algorithm

The physics-based algorithm is founded on the equation of motion of rigid bodies as well as other auxiliary models describing collision phenomena. This algorithm can be described in four main stages: preprocessing, initialization, simulation, and termination.

In the preprocessing stage, the details of rigid objects, i.e., catalyst pellets and container, including the physio-mechanical properties as well as their 3d models in the world space, are defined. This (preset) information is summarized in Table .

Table 1

Preset Data in Preprocessing Stage

The initialization stage involves the initial placement of the pellets at the top of the tube as well as setting all thresholds for friction and collision/contact detection subalgorithms. Having placed the pellets at the top of the container, we set the time interval Δt for the iterative calculation procedure (in our tests 1/40 s) as well as the desired total simulation time. We can then proceed with the simulation, in which the ODE solver is run simultaneously with the collision search subalgorithms to determine the trajectories of the barycenter and orientations of the pellets in the bed at each time step. To inspect whether dynamic equilibrium occurs, i.e., the termination stage, the work-energy theorem is exercised: a system of moving pellets reaches a dynamic equilibrium if the change in total kinetic energy of the system over a typical time period approaches zero. Figure exhibits the skeleton of the proposed packing algorithm.

Figure 6

Flowchart of the packing algorithm.

Results and Discussion

Several series of packing simulations were conducted to assess the influence of preset parameters, i.e., physio-mechanical properties such as coefficient of restitution (COR) and surface friction factor, as well as initial placement of pellets at the top of the tube, on the structural properties of the packings. We investigated packing structures of spheres, solid cylinders, and Raschig rings, with a tube-to-pellet diameter ratio, N, within the range 3.06–9.16.

To inspect the fidelity and robustness of the proposed packing algorithm, the simulated packings have been scrutinized and benchmarked, first and foremost, in terms of bulk porosity and radial void fraction distribution. To this end, the experimental data from refs (25), (28), and (76) for bulk porosity and refs (24), (26), (27), (52), (77), and (78) for radial void fraction distribution were tested, which all together provide data for a wide range of tube-to-pellet diameter ratios, 1.7 ≤ N ≤ 27.95.

Roles of Loading Methods and Filling Speed

The choice of initial placement of catalyst pellets, i.e., the loading method, in the setup of the packing simulations is of paramount importance as it may have a considerable influence on the dynamic behavior of the packing process and densification, and accordingly, on the topological properties of the generated structures. The packing algorithm is set up on the basis of the data given in Table , with COR = 0.6 and μd = 0.6, to simulate random packings of spheres with three different loading strategies: (1) in the first scheme the pellets are introduced in a column in line with the tube axis; (2) in the second scheme the pellets are placed in two colums, equidistant from the tube axis; (3) in the third scheme the pellets are placed in four columns, equidistant from the tube axis. These scenarios (see Figure ) have been examined for two packings with tube-to-pellet diameter ratios of N = 3.1 and 6.1.

Figure 7

Various scenarios considered for loading method: (a) spherical packing with N = 3.1; (b) spherical packing with N = 6.1.

The simulated structures then were examined in terms of bulk voidage to determine which scenario would result in a denser structure. The mean voidage of the generated packings was computed precisely on the basis of the mesh counts of catalyst pellets up to a packing altitude of H = 100 mm and were benchmarked against empirical correlations proposed by Dixon[25] and Foumeny et al.[76] The results are shown in Figure .

Figure 8

Influence of different loading schemes on the bulk voidage of generated structures: (a) spherical packing with N = 3.1; (b) spherical packing with N = 6.1.

As demonstrated in Figure , the first loading strategy not only results in denser packings among all test cases but also is in better agreement with empirical correlations describing dense packing structuers. This can be attributed to the higher intensity of lateral displacements of the pellets for loading method 1, leading to a higher order of densification. The same analysis has been performed for packings of cylindrical pellets with N = 6.98; the results of the mean porosity analysis have confirmed the superiority of the first scenario, where the cylinders are placed obliquely with an angle of 45° with respect to the gravity direction. This leads to the densest structure for such a narrow fixed beds and is in the best agreement with published correlations. Similar results have been found by Fernengel et al.[79] for random packing of spheres with N = 6.25, where the authors investigated the influence of number of spheres per layer in the loading scheme on bulk voidage of simulated structures using LIGGGHTS and Blender packages. It was also found that increasing the distance between successive pellets in their initial placement allows the tube to be a little more compacted. In a second step, we imposed random translational and orientational disturbances on the pellets in the first scheme of loading, resulting in five extra loading schemes, to account for the influence of random pouring on the bulk voidage of RBD-simulated structures. This analysis was conducted for random packings of spheres (including translational disturbances) with N = 3.1 and 6.1 and cylinders (including both translational and orientational disturbances) with N = 3.55 and 6.98. Figure illustrates the imposed disturbances on the first loading scenario for packing of cylinders with N = 6.98.

Figure 9

Random disturbances imposed on the first loading scheme to mimic the influence of random pouring on bulk voidage in random packings of cylinders with N = 6.98.

Overall, the results show that denser structures can be synthesized using the original version of the first loading scheme for all cases. Here, for sake of brevity, the result of comparison analysis for random packing of cylinders with N = 6.98 is presented in Table , where the bulk voidage obtained from the generated packings according to the loading schemes presented in Figure is compared with the empirical coloration by Dixon.[25]

Table 2

Influence of Imposing Translational and Orientational Disturbances on the First Loading Scheme on the Resulting Bulk Voidage for Packing of Cylinders with N = 6.98

	loading scenario
	1	2	3	4	5	6
bulk voidage (calculated)	0.405	0.410	0.419	0.416	0.420	0.425
predicted voidage after Dixon (1988)	0.395	0.395	0.395	0.395	0.395	0.395
MER (%)	–2.4	–3.8	–6.1	–5.2	–6.2	–7.5

It is worth remarking that the aforementioned loading schemes were also implemented in our further analysis described in sections and 5.3, where we have investigated the roles of restitution and friction factors in the bulk voidage of simulated structures. Similarly, our results demonstrate that denser packing structure can be generally generated using the original version of first loading method for such narrow beds for all ranges of restitution and friction coefficients considered.

Role of Surface Bounciness (Coefficient of Restitution)

The amount of surface bounciness, expressed as the coefficient of restitution, can substantially affect the subsequent chain of collisions, and therefore the final structures. Zhang et al.[80] investigated the influence of the coefficient of restitution on packing density, but their work was restricted to packings of spheres in large-N beds, e.g., N = 24, where in the effect of container wall on the packing process and final bulk voidage cannot be elucidated. Since the impact of this physio-mechanical property on the packing densification has not yet been systematically investigated in low-N packings of (non-) spherical pellets, an effort is made here to monitor the role of this property on the bulk porosity of packings of spheres and cylinders. To this end, the packing algorithm is set up on the basis of the data given in Table , to synthesize random packings of spheres with N = 3.1 and 6.1 and cylinders with N = 3.55 and 6.98. The friction factor is set to μd = 0.6, while the restitution coefficient is set to 0.05, 0.15, 0.35, 0.55, 0.75, and 0.95, respectively, covering the whole range of collision behavior from semiplastic to semielastic. The first loading scheme is applied. However, to generate statistically independent samples for averaging, random disturbances (with a maximum magnitude of 0.8 dp and maximum angle of 45°) have been imposed to the initial pellet positions and orientations for each packing case (see Figure ). This was repeated six times for each case, leading to 144 generated test cases for spherical and cylindrical packings.

The mean voidage of all test cases was then computed on the basis of the number of pellets stacked within a tube up to the altitude of H = 120 mm and depicted versus COR in Figure .

Figure 10

Influence of coefficient of restitution (COR) of catalyst pellets on the bulk porosity of generated structures for packings of (a) spheres with N = 3.1, (b) cylinders with N = 3.55, (c) spheres with N = 6.1, and (d) cylinders with N = 6.98. Results of six independent simulations for each COR (with different random disturbances on the initial pellet positions and orientations) are shown, indicating the variability of the results. Dashed lines are trendlines based on the average of these independent simulations.

Overall, the results of this analysis agree with the prevailing hypothesis on the influence of restitution coefficient, viz., higher values of COR generally result in denser packings, because with more elastic collisions, the probability of longer-lasting successive collision chains increases. This allows pellets to be further displaced laterally, and to be vibrated for a longer duration of time, allowing them to find their optimal positions. However, our results also show considerable variability of the mean porosity from sample to sample, in particular in narrower structures (see Figure a,b), which can be reasonably attributed to the restrictive role of the tube wall. The trends in wider beds, as shown in Figure c,d, demonstrate a lower amount of variability. Furthermore, the results demonstrate that the influence of COR on the resulting packing density is more discernible for COR values beyond 0.5, where it causes an intensive vibration of catalyst pellets in the bed.

Role of Surface Roughness

The surface roughness of a pellet, which is fundamentally described by its friction coefficient, is another physio-mechanical property that can affect the process of packing. To assess the impact of this parameter on structural properties of random packings, the same procedure as for the role of COR is pursued. In this case, the pellet’s COR is set to 0.6, whereas the surface friction coefficient is set to 0.05, 0.3, 0.5, 0.7, and 0.95, respectively, covering a large spectrum of dynamic friction coefficient.

The bulk porosity of all samples for each test case has been computed, and plotted against friction coefficient in Figure .

Figure 11

Influence of friction factor of catalyst pellets on the bulk porosity of generated structures for packings of (a) spheres with N = 3.1, (b) cylinders with N = 3.55, (c) spheres with N = 6.1, and (d) cylinders with N = 6.98.

The results confirm the general understanding of the effect of surface friction on packing density: lower values of friction facilitate sliding of the contact surfaces relative to each other, leading to denser packing structures. In fact, the influence of friction factor on the packing density appears to be 1.5 times larger (in the studied range) than that of the COR, leading to an apparently lower variability between different samples.

Validation Study and Postprocessing of the Results

The main goal of our work is to introduce a new RBD methodology to synthesize realistic random packings of particles of any shape. We now present a detailed validation study, whereby the reliability of the RBD-algorithm in replicating structural properties of realistic packings is scrutinized. The specifications of the physical simulations (see Table ) were chosen to reproduce the experimental arrangements utilized in the works of Benenati and Brosilow,[24] Roshani,[77] and Mueller.[52] More simulation runs have been performed to generate random packings in the range 3.1 ≤ N < 9.16 to examine the restrictive role of confining walls on the predicted bulk porosity over a range of tube-to-pellet diameter ratios, and to benchmark the predicted radial porosity distribution against some of the most-frequently used empirical correlations. Table addresses the specifications and data used in setup of RBD simulations.

Table 3

RBD Simulation Setup for Validation Study

catalyst pellet	bed environment
no. of face mesh (per pellet):	bed size:
sphere: 3120	spherical beds: dt = 31, 39.6, 41, 56, 59.6, 61, 79.9 mm & bed altitude: 120 mm
cylinder: 4400	cylindrical beds: dt = 33.65, 35.49, 46.93, 69.83, and 91.58 mm & bed altitude: 120 mm
Raschig ring: 8008	Raschig ring beds: dt = 30.58, 40.45, 60.18, and 79.91 mm & bed altitude: 120 mm
density: 8030 kg/m3	tube wall friction coefficient (dynamic): 0.6
surface friction coefficient (dynamic): 0.1	tube wall surface bounciness (COR): 0.6
surface bounciness (COR): 0.9	gravity acceleration: 9.81 ms–2

Figure illustrates the dynamic behavior of the random packing process, in four frames, for a typical cylindrical packing. The Supporting Information contains a video of a similar packing process.

Figure 12

Dynamics of packing process for packing of cylinders with COR = 0.9 and μd = 0.1 in a tube with N = 4.69. A video of a packing process with COR = 0.1 and μd = 0.9 with N = 5 is supplied in the Supporting Information.

Several packings have been generated, on the basis of the specifications in the Table , including seven different tube-to-particle ratios for packings of spheres, five for packings of solid cylinders, and four for packings of raching rings. Each case has been executed in triplo with slightly perturbed initial pellet positions to generate statistically independent data. Figure illustrates some typical results of RBD-simulated packings of spheres, cylinders, and Raschig rings.

Figure 13

RBD-simulated structures.

The bulk porosity of the generated structures has been computed on the basis of the total amount of pellet material up to the altitude of 100 mm and compared with well-known published correlations in Figure .

Figure 14

Comparison between the mean porosity extracted from RBD simulations and empirical correlations for three types of pellets: (a) spheres, (b) equilateral solid cylinders, and (c) Raschig rings.

The RBD simulation results demonstrate satisfactory agreement with the empirical correlations for all cases, giving a maximum relative error (MRE) of 3.9%, 17.6% (relating to the narrowest bed case), and 6.6% for packings of spheres, cylinders, and Raschig rings, respectively, based on the empirical correlation by Dixon.[25]

The local structural properties of the RBD-simulated packings have also been examined. For this, the axially averaged radial void fraction distribution of the RBD-simulated structures was extracted and compered to literature data. To evaluate radial void fraction profiles, a planar mesh-based approach was adopted, in which a packing structure is intersected with a series of concentric tubes with different diameters (see Figure ). The local void fraction at a specific radius can then be precisely computed by the following formula:where Sint is the intersecting area with pellets and Stotal is the area of tube that intersects the packing at radius r.

Figure 15

Implementation of the planar mesh-based approach for evaluating axially averaged radial void fraction data for RBD-simulated structure of raching ring with N = 6.02.

The radial void fraction profiles from the RBD simulations versus the distance from the tube wall (made dimensionless by dpv, the equivalent diameter of a sphere of the same volume) are shown in Figures and 17 for spheres and cylinders, respectively.

Figure 16

Comparison between the radial void fraction profiles obtained from RBD simulations of sphere packings and literature data. The experimental data presented here are extracted with permission from refs (52) and (24).

Figure 17

Comparison between the radial void fraction data obtained from RBD simulations of packings of cylinders and experimental data by Roshani[77] and Giese et al.[78] The experimental data presented here are extracted with permission from refs (77) and (78).

Figures and 17 show very good agreement between the radial void fraction profiles extracted from RBD simulations and the corresponding experimental data by Mueller[52] for packings of spheres and Roshani[77] and Giese et al.[78] for packings of cylinders. The same conclusion can be made when the RBD data are compared with empirical correlations by de Klerk[27] and Roshani,[77] giving R2 and root-mean-square error values of more than 0.85 and lower than 0.08 for spheres, respectively, and more than 0.84 and lower than 0.07 for cylinders, respectively.

Overall, it can be concluded that the proposed packing algorithm is able to reasonably reproduce the essential features, including the oscillatory-damped behavior as well as amplitude and period of oscillations, of the most-frequently used experimental data.

We can now use the method to investigate the radial void fraction profile in packings of Raschig rings, as shown in Figure .

Figure 18

Radial void fraction profiles obtained from RBD simulations of packings of Raschig rings.

As revealed in Figure , the behavior and pattern of the void fraction distribution over the dimensionless distance from the tube wall in Raschig ring packings is quite different compared to those known for sphere and cylinder packings. A similar trend has been reported by Giese et al.[78] for a Raschig ring bed with N = 10.

Conclusion

A physics-based packing algorithm, founded on the concepts of rigid body dynamics (RBD), has been presented and validated. The algorithm enables us to synthesize random packing structures of nonspherical and nonconvex shapes. The advantage of our approach, compared to popular computer graphics software such as Blender, is that we proposed a more rigorous and realistic approach to model the resting contact phenomenon, wherein the transition between moving and resting particles is controlled by a cutoff on the relative contact velocity followed by a detailed balance of constrained forces acting on each pellet, facilitating the stability of convergence in RBD simulations. Our approach avoids the usage of artificial translational and angular momentum sinks on moving particles and therefore exactly obeys the laws of conservation of linear and angular momentum.

We used our new approach to generate realistic random packings of spheres, cylinders, and Raschig rings with a bed-to-pellet diameter ratio N ranging from 3 to 9.16, where the role of confining walls in the packing process is very important. The results of our validation study have demonstrated satisfactory agreement with literature data concerning the bulk porosity and radial void fraction profiles, substantiating the merits of this approach in replicating the structural properties of realistic random packing geometries of nonspherical catalyst pellets. Furthermore, the influence of essential physio-mechanical properties of catalytic particles, such as surface roughness and bounciness, were studied for both spherical and cylindrical particulate beds.

A particular feature of this approach is that it can synthesize both loose and dense packing structures, depending on the precise loading scheme and physio-mechanical properties, thereby in a way mimicking both sock-loading and dense-loading methods utilized in industrial practice.

The packing algorithm provides detailed information concerning the topological features of randomly packed fixed bed structures, e.g., the position and orientation of catalyst pellets in the bed. Therefore, it has the potential of being used as a supplementary tool to lattice Boltzmann or computational fluid dynamics simulations of reacting flows and heat- and species-transport characteristics of such complicated unit operations.

We finish our conclusions by noting that from the results presented in this work, we cannot yet make a statement about the possibly improved quality of the packing structures by our method relative to other methods such as Blender, DigiDEM, and LIGGGHTS. We have shown that the packing structures are highly sensitive to pellet shape, physiomechanical properties, and loading methods. To trustfully show the differences between the different methods, it is therefore necessary to conduct a systematic comparison study, covering different pellet shapes with different ranges of physiomechanical properties and loading methods. We note that small but significant differences in the void fraction of sphere packings in cylindrical beds predicted by Blender, compared to DigiDEM and LIGGGHTS, have already been found in a recent paper by Fernengel et al.[79] (see their Figure 2). At this point it is not certain whether such differences are caused by differences in handling the dynamics of the packing process (which for instance is influenced by the magnitudes of the translational and angular momentum damping terms) or by small differences in the contact model. Small differences in predicted voidage may even be acceptable for certain applications in view of the greatly enhanced computational speed of Blender relative to these other packages. In any case, a careful and systematic comparison of these different methods is clearly needed and will be the topic of our forthcoming contribution.