**1. Introduction**

Fixed-bed reactors are heavily used in the chemical and process industry, especially in the field of heterogeneous catalysis, where there are thousands of individual catalytic fixed-bed reactors with a low tube-to-particle diameter ratio *N* ≤ 10 that are interconnected to tube-bundle reactors. This design decision is the result of optimizing multiple objectives, such as low pressure drop, good radial heat transport, and high active catalytic surface area [1]. Nevertheless, advancing climate change and the shortage of raw materials make more resource- and energy-efficient processes necessary. Here, numerical methods can play a paramount role to develop better designs faster and that are more cost efficient. In the last few years, particle-resolved computational fluid dynamics (CFD) was heavily used by numerous authors to develop process intensification strategies with the focus on the effects on the mesoscopic pellet scale. The range of works extends from investigations of the influence of particle shape on bed morphology and fluid dynamics [2–4], heat transport [5–8], and mass transfer processes [9–12] to the development of novel reactor concepts, such as packed foams [13–15], periodic open-cell structures [16–18], finned reactors [19], or the use of random macroscopic wall structures [20,21]. However, particleresolved CFD is a numerically very demanding method, and its applicability is currently limited to systems with a few thousand particles [22].

When talking about process intensification, however, this can take place at different spatial scales [23], ranging from the molecular scale to the pellet scale to plant scale. On

**Citation:** Jurtz, N.; Srivastava, U.; Moghaddam, A.A.; Kraume, M. Particle-Resolved Computational Fluid Dynamics as the Basis for Thermal Process Intensification of Fixed-Bed Reactors on Multiple Scales. *Energies* **2021**, *14*, 2913. https://doi.org/10.3390/en14102913

Academic Editors: Goodarz Ahmadi, Kiao Inthavong and Pouyan Talebizadeh Sardari

Received: 31 March 2021 Accepted: 11 May 2021 Published: 18 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the plant-scale level, process intensification options are process optimization [24–26], the development of process integration concepts [1,27,28], and applying dynamic operating conditions [29,30]. Due to the restrictions discussed above, particle-resolved CFD cannot directly be applied for the simulation on the largest scale. For this, process simulation software, e.g., *Aspen Plus*, *gPROMS*, or the open-source solution *DWSIM*, is the more efficient choice. Most often, these software packages use two-dimensional pseudo-homogeneous models to describe fluid dynamics and the heat and mass transfer of fixed-beds. For these kinds of models, the knowledge of effective transport parameters, e.g., the effective viscosity and thermal conductivity, the wall heat transfer coefficient, and the axial dispersion coefficient, is necessary to obtain accurate results. Since those parameters need to lump a series of effects that have their fundamentals on micro- and meso-scale fluid dynamic effects, published data varies greatly [31] and can often only be found for certain reactor designs. This is where particle-resolved CFD comes into play, since it has the potential to act as a data source to derive the effective transport parameters that are needed for a reliable process simulation. Recent publications by Dixon [32] and Moghaddam [33] show encouraging results.

In the scope of this work, we investigated different fixed-bed reactor concepts numerically, using particle-resolved CFD. Besides reactors filled with different particle shapes, namely spheres, cylinders, Raschig rings, and four-hole cylinders, additionally, the impact of macroscopic wall structures was studied for packings of spherical particles. The research focus lied on the quantification and qualitative characterization of their heat transport characteristics. For the sake of reduced complexity and to nail the investigations down to the impact of fluid dynamic effects only, no chemical reactions were considered in the investigated cases. The aim of this study was to:


#### **2. Materials and Methods**

In this section, the fundamentals of the numerical models are briefly discussed. For a more detailed description, the reader is referred to literature that explains the fundamentals of CFD [34], particle-resolved CFD [22] and simplified fixed-bed reactor modeling [24,35] in more detail.

#### *2.1. Particle-Resolved CFD*

Particle-resolved CFD is an established numerical method for the simulation of fixed-bed reactors. This CFD-based modeling approach is characterized by a full threedimensional spatial resolution of all particles and their interstices. The general procedure consists of four fundamental steps: packing generation, CAD generation, meshing, and the CFD simulation itself. For a more detailed discussion about particle-resolved CFD, the interested reader is referred to comprehensive review articles by Dixon et al. [36,37] and Jurtz et al. [22]. A description of all numerical methods used in the scope of this work can be found in Supplementary Material, Sections S1 and S2 attached to the article. The most important material properties and boundary conditions are summarized in Table 1. All numerical simulations were conducted with the commercial CFD tool *Simcenter STAR-CCM+* provided by *Siemens PLM Software*.


**Table 1.** Material properties and boundary conditions for the DEM and CFD simulations.

#### 2.1.1. Numerical Packing Generation Using DEM

In this study, the discrete element method (DEM) was used to numerically generate random packings of spherical and various cylinder-like particle shapes. For the nonspherical particles, the contact detection algorithm of Feng et al. [39] for cylindrical particles was used. The particle beds of Raschig rings and four-hole cylinders are identical to the ones of the cylinders in terms of particle position and orientation, since they are based on the same DEM simulation results. This means that for particles with inner voids, as Raschig rings and multi-hole cylinders, the effect of the inner voids on the particle dynamics during the filling process was neglected. In our previous studies, it was shown that this was a valid assumption [4]. For all filling simulations, the linear spring contact model was used. An overview of all generated packings is given Figure 1.

The aim of this study was to analyze the impact of particle shape and packing mode on the fluid dynamics and heat transfer. Therefore, fixed-beds filled with different particle shapes were generated, whereby the tube-to-particle diameter ratio was held fixed to a value of *N* = 5 by setting a constant volumetric sphere-equivalent particle diameter of *d*p,v = 11 mm and inner tube diameter of *D* = 55 mm. It is known from previous studies [40,41] that even or odd numbered tube-to-particle diameter ratios can lead to additional heterogeneities in the bed morphology. It was expected that additional morphological heterogeneities would lead to an increase of thermal heterogeneities as well. In order to identify the limitations of the pseudo-homogeneous two-dimensional plug flow model, we decided to benchmark the model for such an extreme case against particle-resolved CFD results. Increasing particles thermal conductivity to extreme values, as recently done by Moghaddam et al. [33], is also an option to increase thermal heterogeneities. However, since in most applications, the ceramic-type catalyst support, which is often porous, is used, which is characterized by a low thermal conductivity (*λ*<sup>s</sup> ≈ 0.2W/(m K)), we decided to not choose that option. A sketch of each investigated particle shape, including its dimensions, is given in Figure 2. The bed height was set to *h* = 100 *d*p,v. For each particle shape, two different packing modes were created: a rather loose and a dense bed configuration. To achieve different packing densities, the static friction coefficient was used as a tuning parameter during the DEM simulation. A more detailed description of this method can be found in our previous publications [4,42].

The macroscopic wall structure was generated with a Java macro. Along the reactor length of 1.1 m, on 96 axial stages, fifteen spheres were homogeneously distributed on each stage, whereby the center of mass of each sphere was in coincidence with the tube radius. Subsequently, each sphere was move in the outward direction by a factor of rand(0, 1) · *d*p,v/2. Thereafter, all spheres of each stage were rotated with a random angle around the reactor axis, using a rotation angle of rand(0, 1) · 2*π*. In a final step, the macroscopic wall structure was generated by subtracting all spheres from the reactor tube. The structured tube was afterwards filled with spherical particles, whereby the same simulation parameters and boundary conditions were used as for the smooth wall setup.

**Figure 1.** Overview of all generated packings. Left: loose packings of (**A**) spheres, (**C**) cylinders, (**E**) rings, (**G**) 4-hole cylinders, and (**I**) spheres with a macroscopic wall structure. Right: dense packings of (**B**) spheres, (**D**) cylinders, (**F**) rings, (**H**) 4-hole cylinders, and (**J**) spheres with a macroscopic wall structure.

**Figure 2.** Investigated particle shapes: (**A**) spheres, (**B**) cylinders, (**C**) rings, and (**D**) 4-hole cylinders.

#### 2.1.2. Meshing

After the bed generation was completed, the position and orientation of all particles were extracted, and based on these data, a CAD model of the fixed-bed was generated by placing CAD parts of the respective particle shape. Subsequently, the geometry was meshed, whereby the improved local "caps" approach, developed by Eppinger et al. [8], was used to avoid bad cell quality near particle–particle and particle–wall contacts. This

enhanced meshing strategy was based on earlier work of the authors [41]. In one of our previous studies, it was proven that this meshing approach not only worked fine for spherical particles, but also for more complex particle shapes such as cylinders, Raschig rings, and multi-holed cylinders [4,43]. It was found that the mesh settings used led to mesh-independent results regarding fluid dynamics and heat transfer [3,10,44]. Recently, Eppinger and Wehinger [8] investigated the impact of the gaps that were introduced between the particles during the meshing process. They found that the gap size had only a marginal effect on bed voidage and pressure drop. However, the authors found that the fluid in the gaps was no longer stagnant if the gap size was increased above a value of 0.01 *d*p,v, which was the value that was used in the present work. Therefore, a bigger gap size than the one used in this study could negatively affect the accuracy of heat transfer simulations, since an additional thermal resistance would be introduced for the inter-particle heat transfer.

To avoid unwanted inlet and outlet effects, the inlet and outlet faces were extruded a distance of 1 *D* and 3 *D* away from the bed, respectively. Two prism layers with a target thickness of 0.025 *d*p,v were used to capture the fluid dynamic and thermal boundary layer at the particles and the tube wall.

#### 2.1.3. CFD Simulation

The momentum, energy, and turbulence transport equations were solved in a segregated manner (see Section S2). For the sake of reduced complexity and to be able to study the convective and conductive heat transfer without any superposing heat transfer mechanism, the heat transfer simulations were conducted under conditions where radiative heat transfer can be neglected, which was, therefore, not accounted for. Because of this, an inlet temperature of *T*<sup>0</sup> = 20 °C and a constant wall temperature of *T*<sup>w</sup> = 200 °C were used along the fixed-bed region.

For all simulations, the SIMPLE-algorithm was used for pressure-velocity coupling, and turbulence was considered through Reynolds-averaged Navier–Stokes (RANS) equations in conjunction with realizable *k*-*ε* model-based closures. This turbulence model was successfully used in our previous works [4,8,10,41,43].

#### *2.2. Simplified Heat Transfer Modeling*

The class of pseudo-homogeneous models is widely used for the simulation of fixedbed reactors. Here, the particle scale is not resolved. Instead, all effects are lumped into effective transport parameters. In terms of heat transport, the effective transport parameters needed are either a radially invariant effective thermal conductivity *λ*eff,r and a wall heat transfer coefficient *α*<sup>w</sup> or only a radially varying effective radial thermal conductivity *λ*eff,r(*r*). The two different concepts are described in the following sections.
