Validation of pressure drop prediction and bed generation of fixed-beds with complex particle shapes using discrete element method and computational fluid dynamics

Catalytic fixed-bed reactors with a low tube-to-particle diameter ratio are widely used in industrial applications. The heterogeneous packing morphology in this reactor type causes local flow phenomena that significantly affect the reactor performance. Particle-resolved computational fluid dynamics has become a predictive numerical method to analyze the flow, temperature, and species field, as well as local reaction rates spatially and may, therefore, be used as a design tool to develop new improved catalyst shapes. Most validation studies which have been presented in the past were limited to simple particle shapes. More complex catalyst shapes are supposed to increase the reactor performance. A workflow for the simulation of fixed-bed reactors filled with various industrially relevant complex particle shapes is presented and validated against experimental data in terms of bed voidage and pressure drop. Industrially relevant loading strategies are numerically replicated and their impact on particle orientation and bed voidage is investigated.


| INTRODUCTION
Syngas is feedstock for many important chemical processes like methanol and ammonia synthesis, the production of aldehydes by oxo synthesis, and the generation of hydrocarbons by the Fischer-Tropsch process. Common routes for syngas production are heterogeneous catalytic reforming processes, for example, steam and dry reforming of methane. 1 Because of the highly endothermic or exothermic nature of these reactions, tubular reactors consisting of numerous fixed-beds with a low tube-to-particle diameter ratio N are widely used. They are characterized by an intensified radial heat transfer and a low pressure drop. However, in this reactor type the assumption of a homogeneous radial void fraction distribution is not valid and local flow phenomena play an important role. This leads to a strong interplay between bed morphology, fluid dynamics, heat and mass transfer, and therefore, the reactor performance itself.
When it comes to process intensification, the probably most obvious way to improve the overall process is the testing of different catalyst shapes. The ideal particle shape has to satisfy multiple objectives: a high active catalytic surface, low pressure drop, low axial dispersion, and a good radial heat transfer characteristic. Several authors investigated the effect of the particle shape on reactor performance and life expectation theoretically and experimentally. Bruno et al 2 compared multi-hole cylinders with Raschig rings theoretically and found for multi-hole cylinders a two times longer maintenance interval until the catalyst needs to be replaced and a significantly increased life expectation for the reactor tubes. In their comprehensive review paper Sie and Krishna 3 discuss several aspects regarding the impact of particle shape on reactor performance and found that the right choice of particle shape is a multi-objective optimization problem to solve the complex trade-off between low pressure drop, high surface-to-volume ratio, high specific reaction rate, low manufacturing cost, and high particle strength. Afandizadeh and Foumeny 4 give in their work special attention to the impact of particle shape on bed morphology. They compare spheres, cylinders, Raschig rings, Lessing rings, and cross web cylinders regarding their active surface per unit volume while particle aspect ratio and inner-to-outer diameter ratio is varied. For cylinders they suggest an aspect ratio of (h/d p ) = 1, while for Raschig rings, Lessing rings, and cross web cylinders an optimum design range slightly reducing the conversion of methane steam reforming in comparison to 7-hole cylinders. The latter particle shape used for steam reforming was also investigated by Franczyk et al. 7 They studied the influence of aspect ratio, hole diameter and tube-to-particle diameter ratio on relative activity, relative gas load, pressure drop, and wall temperature theoretically. The authors propose an aspect ratio of one, a tube-to-particle diameter ratio N ≥ 5, and a hole-to-cylinder crosssection area ratio of 0.30-0.37. Karthik and Buwa 8 used particleresolved computational fluid dynamics (CFD) to investigate the impact of different particle shapes on the reactor performance for four industrially important solid-catalyzed gas-phase reactions. They found that, with increasing particle surface area, pressure drop increases while intra-particle temperature and concentration gradients decrease. The authors show that particle shape significantly impacts reactor performance. However, the overall efficiency depends on several factors as mass transfer limitation and reaction equilibrium. Therefore, the optimal particle shape was found to be dependent on the chemical reaction that is considered.
The relatively low number of experimental work in that field is explainable by the high costs and required time that is needed to conduct the experiments. Furthermore, the knowledge gain by experimental results is limited, since they are most often based on inletoutlet measurements only. Therefore, in the last twenty years many researchers developed and used particle-resolved CFD simulations to get a better understanding of fluid dynamics and transport processes inside fixed-bed reactors. Particle-resolved CFD is an almost first-principle modeling approach where the whole reactor including the interstices and if needed also the particles are spatially discretized and conservation equations are solved. Comprehensive reviews of this method can be found in Dixon et al 9 and Jurtz et al. 10 For simple particle shapes like spheres, cylinders, and Raschig rings many authors validated the particle-resolved CFD method against experimental results or correlations in terms of bed voidage, [11][12][13][14][15] pressure drop, 14,[16][17][18] and temperature and species profiles. [19][20][21] However, for more complex particle shapes only few publications exists. Caulkin et al generated packings of 4-hole cylinders, 22 pall rings, 13 and trilobes 23 numerically and compared the radial void fraction distribution with experimental data. For 4-hole cylinders and trilobes the numerical results were in good agreement with experimental data. For pall rings the void fraction was underestimated and only a fair agreement was found. Boccardo et al 16  In order to make particle-resolved CFD a design tool for new particle shapes the currently existing lack of experimental validation needs to be fixed. In this work, packings of different complex particle shapes with inner voids and outer surface structure are numerically generated using DEM and the pressure drop is calculated by CFD simulations. It is shown that a precise representation of the bed morphology is basis for a predictive CFD simulation. Therefore, special emphasis is given to the aspect of numerical packing generation and a method is presented that improves the match between experimentally and numerically generated packings by calibrating the static friction coefficient. 24 In industrial applications, very often special loading devices are used when the reactor is filled with particles. For the first time the impact of a special loading device, which is often used in industry, is considered and its impact on particle orientation and bed voidage is investigated.

| Experimental measurements
Experimental measurements were conducted in cylindrical reactor tubes with an inner diameter of D = 101.6 mm and D = 152.4 mm. Four different particle shapes were investigated. The shapes and their dimensions are depicted in Figure 1 and Table 1. Besides Raschig rings (RR) with almost equilateral dimensions, three more complex shapes were studied: two different types of 10-hole cylinders low pressure drop particle A loading procedure that is often applied in industrial applications was used: as depicted in Figure 2a a CATCADE™ reformer loading device by Cat Tech ® is used for filling. The device is placed in central position inside the reactor tube, as shown in Figure 2b, and particles are poured from the top into the reactor until half of the tube is filled. Afterwards, the bed is densified by hammering against the reactor wall for 30 s as shown in Figure 2c. This procedure is repeated once. During the loading process the CATCADE™ device is moved upwards to prevent that it immerses into the bed.
The loading device is used to reduce particle velocities and to prevent them from damage during the filling process. Furthermore, it is designed to reduce the variance in packing morphology, and therefore, variance in pressure drop between different reformer tubes. The overall aim is a reduced turnaround time for the filling process in industrial applications. The bed voidage is determined by counting the particles that are needed to reach a specified bed height: measuring the pressure difference before and after an orifice plate using HD750 differential pressure manometer by Extech Instruments.
The pressure drop across the sample bed is measured using HHP-2081 digital manometer by Omega Engineering with a maximum uncertainty of ≈30 Pa. The pressure drop was measured between inlet and front of a perforated plate, which is where the particles rest on. For each particle shape and inlet velocity the pressure drop measurements were repeated two times, whereby the bed was newly generated each time.

| Numerical methods
In the past years, several workflows have been proposed by different authors to conduct particle-resolved CFD simulations of fixed-bed reactors, as recently reviewed by Jurtz et al. 10 What all workflows have in common is that they are based on the following four sequential steps: packing generation, CAD creation, meshing, and the actual CFD simulations. In this recent study, we use the approach developed by Eppinger et al. 14,15 Here the bed is created by using DEM. Afterwards, position and orientation vectors of all particles are extracted and stored in a csv file. Based on this data, a CAD description of the packed bed is generated by placing CAD parts of the respective particle shape using JAVA macro functionality. As a final step the geometry is meshed. To avoid bad cells near particle-particle and particlewall contacts, the local "caps" strategy is used 14  Simcenter STAR-CCM+ provided by Siemens PLM Software.

| Numerical packing generation
For numerical packing generation the discrete element method (DEM) established by Cundall and Strack 25 is used. This approach allows particles to overlap to a certain degree to calculate restitution and damping forces based on the overlap. The linear momentum equation for each particle is given by Newton's law of motion: Here m p and υ p are particle mass and velocity, t is time and F s and F b are the sum of surface and body forces that act on the particle.
For the filling process only the gravitational force and the contact forces are considered. The overall contact force is the sum of particleparticle and particle-wall contact forces acting on a particle, whereas each contact force can be decomposed in a tangential (F n,i ) and a normal (F t,i ) acting component: The force in normal direction is defined by: Here K n is the normal spring stiffness, d n the overlap in normal direction, N n the normal damping and v n the normal velocity component of the relative sphere surface velocity at the contact point. The force in tangential direction is defined as: where K t is the tangential spring stiffness, d t the overlap in tangential direction, N t the tangential damping, v t the tangential velocity component of the relative sphere surface velocity at the contact point and C fs the static friction coefficient.
In this study the non-linear Hertz-Mindlin contact model is used.
This leads to: where N n,damp and N t,damp are the normal and tangential damping coefficients that are calculated from the normal and tangential restitution coefficient C n,rest and C t,rest as follows: M eq , R eq , E eq , and G eq are the equivalent values of mass, radius, Young's modulus, and shear modulus of particles A and B during the collision process: where υ is the Poisson ratio.
Besides the linear, the angular momentum of each particle is conserved as well by: where I p is the particle moment of inertia, ω p the angular velocity, r p the position vector from particles center of gravity to the contact point and M c the acting moment due to rolling resistance: Here C fr is the rolling friction coefficient.
In its original form DEM is only applicable for the simulation of spherical particles. An approach to overcome this limitation is the use of composite (or glued) particles. By using the composite particle method, particles with complex shapes are approximated by an arrangement of spheres that are rigidly glued together. The composite particle representations of the particle shapes used in this study can be seen in Figure 1. A number of one-hundred spherical particles is used to approximate the real catalyst shape.
The particles inner voids are neglected for the filling simulation since it can be assumed that their impact on the particle dynamics is low. This is justified because of the small inner hole diameter which ensures that other particles cannot slide into them, and the axisymmetric nature of the investigated particle shape. Therefore, the center of mass of the pellet is not changed by neglecting the inner holes.
Nevertheless, the particle mass is slightly increased if the inner holes are neglected. While particle sedimentation velocity is unaffected by particle mass, it has an impact on the calculated restitution forces during collision events. While the impact of frictional forces on bed voidage is widely accepted, only little data exists on the impact of particle mass or its density. Pottbäcker and Hinrichsen 26 studied experimentally the impact of different material properties on bed voidage and suspected that the coefficient of restitution and the particle density might affect the bed voidage. However, not all data fit to their hypothesis, for example, for spheres of rusty steel and polyoxymethylene the resulting bed voidage is equal, although, particle density is more than five times higher for steel, and static friction coefficients are similar. This indicates that the impact of particle density on bed voidage might be lower than expected. Hoffmann and Finkers 27 consolidated experimental data and derived a correlation for bed voidage as a function of sphericity, particle diameter, and particle density. However, the important frictional forces were not considered, and their data still shows scattering. Therefore, no clear conclusions can be drawn. To understand the impact of particle density on bed voidage two additional DEM filling simulations were conducted. Five hundred cylindrical particles (h = d p = 16.5 mm) are poured into a cylindrical container with a diameter of 101.6 mm. Randomly distributed injection points were placed at a height of 470 mm where particles enter the domain. Only gravitational and particle-particle and particle-wall contact forces were considered. Since in this case ideally cylindrical particles are considered the sophisticated contact detection method of Feng et al 28 in combination with the Linear Spring Contact Model was used. All boundary conditions and material properties except for the particle density are equal. As reference case a particle density of 2,200 kg/m 3 was used. For the second simulation this value was reduced to 1,100 kg/m 3 . This is comparable to simulating particles with 50% inner voids but neglecting the additional voidage in the DEM simulation. The result of this preliminary study is given in Figure 3. For better visibility the highest and lowest particle position  24 it is possible to obtain a correlation for more complex catalyst shapes with a set of DEM simulations and just one filling experiment. To ensure that the correlation only accounts for frictional and vibration effects the loading device will be considered in these simulations.
To include the CATCADE™ loading device in the filling simulation, as depicted in Figure 4a, the overset mesh approach is used. This method often called chimera grid method was first introduced by  Here ρ is mass density of the fluid, v fluid velocity, and T the stress tensor:

| Computational fluid dynamics
where p is pressure, μ dynamic viscosity, I the unit tensor, and D the deformation tensor: Turbulent effects are considered by using the Realizable k-ε turbulence model which has already successfully been used in previous studies. 14,15,18 The steady-state solution is solved by using an incompressible segregated solver approach. The SIMPLE-algorithm is used to solve for the pressure-velocity coupling.
The numerical domain is spatially discretized using polyhedral cells as depicted in Figure 5. In order to avoid insufficient mesh quality at particle contacts, the caps method introduced by Eppinger et al 14 was applied, which creates small gaps between touching particles. The approach is fully automatized and can be applied to any type of particle shape without additional modification. 10 In comparison to other contact modification methods (e.g., bridging particle contacts) it can directly be applied for heat transfer simulations without having to specify an additional thermal resistance at the contacts. 19 A detailed mesh refinement study was recently published by Minhua et al. 34 Two prism layers are used to resolve velocity gradients at the particle surface and reactor wall. In order to minimize the influence of the boundary conditions at the inlet and outlet, the volume mesh is extended by extruding the inlet and outlet. The total cell count was 9.1 million for Raschig rings and around 11.5 million for LDP and FD shape. The overall mesh quality is reasonably good for such complex geometries. A maximum of one out of ten thousand cells have a least square cell quality below 0.001. This metric is an indicator of the quality of a cell, using the physical location of a cell centroid relative to the cell centroid locations of its face-neighbors. 35 3 | RESULTS

| Pressure drop
The pressure drop of packed beds with a bed height of H = 445.0 mm and a diameter of D = 101.6 mm were investigated numerically and experimentally for Raschig rings, LDP 19 × 16 and FD particle shape. Since pressure drop is very sensitive to bed voidage, it is essential that the difference in bed voidage is as low as possible between the numerical packing and its experimental counterpart. Therefore, the static friction coefficient was adapted using trial-and-error to a value of C fs = 0.02 in the DEM simulation to meet the experimental bed voidage. A comparison between the experimental and numerical values in terms of particle count and bed voidage can be found in Table 2, whereby inner particle voids were considered to calculate the bed voidage. A very good agreement with the experimentally determined particle count can be observed for all three particle shapes. A maximum error of 1.4% can be found for Raschig rings.
In Figure 6 the experimental and numerical packing are visually compared. It can be seen that at the reactor wall the particle orientation is qualitatively in a good agreement. The majority of particles are either aligned orthogonally or in parallel to the reactor wall. There is also a trend of building stacked structures, which is something that Zhang et al 36  and simulation results is given in Figure 7a. The same trend is predicted by experiments and numerical simulations: the largest pressure drop can be observed for Raschig rings followed by the LDP 19 × 16 and the FD particle shape. For the same superficial inlet velocity the pressure drop can be cut in half by using the FD particle shape. The parity plot in Figure 7b shows the accuracy of calculated  to 25%. It cannot be seen that the change in orthogonally aligned particles leads to any different preferred orientation angle. Regardless of whether the loading device is used or not, with increasing friction coefficients the particle orientation distribution gets more uniform, although, particles with an orientation angle below 30 are always underrepresented.
Besides particle orientation also bed voidage is affected by the loading device as it can be seen in Figure 10. limit. With increasing aspect ratio, the effect gets less pronounced due to the increasing amount of kinetic energy that is dissipated by particle-particle friction. Overall, the results show that whenever loading devices are used they should also be considered in the numerical filling simulation since they strongly affect particle orientation and bed voidage.    Figure 11. A bed height of H = 508.0 mm was used for the 6 00 -tube. In case of the 4 00tube the experimental and numerical particle count is almost identical which is not unexpected since this was the reference case for the calibration of the friction coefficient. Nevertheless, the results for the 6 00tube also show very good agreement between experimental particle count and simulation result with an error in terms of particle count below 1%. This indicates that with a once calibrated friction coefficient bed voidage of reactors with different tube-to-particle diameter ratios can be predicted using DEM. This brings particle-resolved CFD one step further toward being a predictive design tool for new catalyst shapes. A potential work flow could be as follows: 1. For a new particle shape a series of DEM simulations with different values for the static friction coefficient are conducted. If needed, a loading device can be included in the DEM simulation to account for its effect explicitly.

2.
A regression curve is fitted based on the numerical data of the bed voidage.
3. A filling experiment with a certain tube-to-particle diameter ratio is conducted and bed voidage is evaluated.
4. Based on the numerical regression curve and the experimentally measured bed voidage the static friction coefficient is calibrated.
5. The calibrated static friction coefficient can be used for filling simulations in the particle-resolved CFD framework to predictively generate packing geometries for different tube-to-particle diameter ratios under same filling conditions.
A good description of bed morphology is of major importance for an accurate prediction of fluid dynamics, heat/mass transfer and chemical reactions by particle-resolved CFD simulations. The proposed work flow significantly reduces the numerical turnaround time and experimental effort to generate packing geometries that are comparable to its experimental/industrial counterpart.
The experiments and simulations were repeated for 19 × 16 LDP shapes using C fs = 0.01 to evaluate if the determined static friction coefficient can also be used for slightly modified particle shapes.
However, the simulation results have shown an up to 9% higher particle count than the one that was found in experiments. Therefore, based on the already existing calibration curve for 19 × 12 LDP given in Figure 10, the static friction coefficient was adapted to a value of  Exp. DEM F I G U R E 1 1 Comparison of experimental and numerical particle count for 19 × 12 LDP (C fs = 0.01) and 19 × 16 LDP shape (C fs = 0.3) for different tube diameters simulations were repeated with the adapted friction coefficient. The comparison to experimental data is given in Figure 11 and shows good agreement between experimentally and numerically generated packings. For the 4 00 -tube the simulation results exactly match the experimental data, while for the 6 00 -tube the particle count is underestimated by 3%. This may indicate that an existing calibration curve, to some extent, can also be used if particle geometry is slightly modified, for example, small changes in aspect ratio. This is in accordance with results of Jurtz et al 24 who found for cylindrical particles only little effect of aspect ratio on bed voidage if 0.5 < h/d p < 2. Nevertheless, further experimental and numerical investigations are needed to assess the validity and accuracy when using a calibration curve for modified particle shapes, especially for those with additional outer surface structure.

| CONCLUSION
Particle-resolved CFD was applied to simulate the fluid dynamics in packed-beds filled with different complex particle shapes over a wide range of industrially relevant Reynolds numbers (Re p ≈ 1300-5800). It was shown that an accurate prediction of pressure drop with a deviation in the range of ±15% can be achieved if bed voidage is in agreement with the experimentally investigated packing. The particle orientation close to the reactor wall was also found to be qualitatively in accordance with lab results.
For the first time the impact of loading devices on particle orientation and bed voidage was studied numerically. It was found that bed voidage is significantly increased if loading devices are used. Therefore, it is advised to explicitly include them in the packing generation simulation. It was shown that the overset mesh approach is an efficient and reliable method to include moving boundaries in the simulation.
Since it is of major importance to meet the bed voidage as close as possible, it was shown how this can be efficiently achieved, and experimental and numerical effort can be minimized. It is proposed to use an adapted static friction coefficient to account for artificial densification and surface characteristics. A once calibrated friction coefficient can be used to simulate reactor configurations that differ in tube-to-particle diameter ratio as long as the same filling strategy is used.