Particle-Resolved Computational Fluid Dynamics as the Basis for Thermal Process Intensification of Fixed-Bed Reactors on Multiple Scales

Abstract

Process intensification of catalytic fixed-bed reactors is of vital interest and can be conducted on different length scales, ranging from the molecular scale to the pellet scale to the plant scale. Particle-resolved computational fluid dynamics (CFD) is used to characterize different reactor designs regarding optimized heat transport characteristics on the pellet scale. Packings of cylinders, Raschig rings, four-hole cylinders, and spheres were investigated regarding their impact on bed morphology, fluid dynamics, and heat transport, whereby for the latter particle shape, the influence of macroscopic wall structures on the radial heat transport was also studied. Key performance indicators such as the global heat transfer coefficient and the specific pressure drop were evaluated to compare the thermal performance of the different designs. For plant-scale intensification, effective transport parameters that are needed for simplified pseudo-homogeneous two-dimensional plug flow models were determined from the CFD results, and the accuracy of the simplified modeling approach was judged.

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\le 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,3,4], heat transport [5,6,7,8], and mass transfer processes [9,10,11,12] to the development of novel reactor concepts, such as packed foams [13,14,15], periodic open-cell structures [16,17,18], finned reactors [19], or the use of random macroscopic wall structures [20,21]. However, particle-resolved 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 the plant-scale level, process intensification options are process optimization [24,25,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:

• understand the effect of particle shape and macroscopic wall structures on the packing morphology and, with this, the fluid dynamics and heat transport in fixed-beds;
• quantify improvements in the fixed-bed reactor design that can be achieved, only from a fluid dynamics point of view;
• increase the phenomenological understanding of fluid dynamics and heat transfer in fixed-bed reactors;
• show how effective transport parameters, such as the effective thermal conductivity and wall heat transfer coefficient, can be extracted from particle-resolved CFD results. Those parameters can then be used in simplified process simulation models for process intensification on the plant scale.

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 three-dimensional 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. Material properties and boundary conditions for the DEM and CFD simulations.

DEM Simulation
Normal/tangential spring stiffness (N/m)	$1\times {10}^6$
Static friction coefficient (-)	0.61/0.01 (loose/dense bed)
Normal/tangential restitution coefficient (-)	0.5
CFD Simulation
Fluid density (kg/m3))	Ideal gas law
Fluid specific heat (J/(kg K))	1006.82
Fluid thermal conductivity (W/(m K))	$0.02414{\left(\frac{T}{273.15\mathrm{K}}\right)}^{3/2}\left(\frac{273.15\mathrm{K}+194.0\mathrm{K}}{T+194.0\mathrm{K}}\right)$
Fluid dynamic viscosity (Pa s)	$1.716\times {10}^{-5}{\left(\frac{T}{273.15\mathrm{K}}\right)}^{3/2}\left(\frac{273.15\mathrm{K}+111.0\mathrm{K}}{T+111.0\mathrm{K}}\right)$
Particle density (kg/m3)	1500
Particle thermal conductivity (W/(m K))	$0.21+1.5\times {10}^{-4}T$ [38]
Particle specific heat (J/(kg K))	1046.7
Inlet velocity (m/s)	0.138; 0.688; 1.376; 2.75
Inlet temperature (K)	293.15
Wall temperature (K)	473.15
Pressure (bar)	1.01325

. All numerical simulations were conducted with the commercial CFD tool Simcenter STAR-CCM+ provided by Siemens PLM Software.

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 non-spherical 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_{\mathrm{p},\mathrm{v}}=11\mathrm{mm}$ and inner tube diameter of $D=55\mathrm{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 (${\lambda }_{\mathrm{s}}\approx 0.2\mathrm{W}/(\mathrm{m}\mathrm{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=100d_{\mathrm{p},\mathrm{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 $\mathrm{rand}(0,1)·d_{\mathrm{p},\mathrm{v}}/2$. Thereafter, all spheres of each stage were rotated with a random angle around the reactor axis, using a rotation angle of $\mathrm{rand}(0,1)·2\pi$. 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.

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.01d_{\mathrm{p},\mathrm{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 $1D$ and $3D$ away from the bed, respectively. Two prism layers with a target thickness of $0.025d_{\mathrm{p},\mathrm{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_0=20°\mathrm{C}$ and a constant wall temperature of $T_{\mathrm{w}}=200°\mathrm{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 fixed-bed 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 ${\lambda }_{\mathrm{eff},\mathrm{r}}$ and a wall heat transfer coefficient ${\alpha }_{\mathrm{w}}$ or only a radially varying effective radial thermal conductivity ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$. The two different concepts are described in the following sections.

2.2.1. Pseudo-Homogeneous ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ Model

The pseudo-homogeneous two-dimensional plug flow heat transfer model under steady-state conditions is described by:

$${\rho }_{\mathrm{f}}{c}_{\mathrm{p},\mathrm{f}}{u}_{\mathrm{z}}\frac{\partial T}{\partial z}={\lambda }_{\mathrm{eff},\mathrm{r}}\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right)+{\lambda }_{\mathrm{eff},\mathrm{z}}\frac{{\partial }^{2}T}{\partial z^2},$$

(1)

whereas the following boundary conditions are used:

$$\begin{array}{ccc}\hfill -{\lambda }_{\mathrm{eff},\mathrm{r}}\frac{\partial T}{\partial r}& ={\alpha }_{\mathrm{w}}\left(T-T_{\mathrm{w}}\right)\hfill & \hfill \mathrm{at}r=R\end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{\partial T}{\partial r}& =0\hfill & \hfill \mathrm{at}r=0\end{array}$$

(3)

$$\begin{array}{ccc}\hfill T& =T_0\hfill & \hfill \mathrm{at}z=0\end{array}$$

(4)

$$\begin{array}{ccc}\hfill \frac{\partial T}{\partial z}& =0\hfill & \hfill \mathrm{at}z=L.\end{array}$$

(5)

Here, $u_{\mathrm{z}}=u_0/\epsilon$ is the constant interstitial velocity, ${\rho }_{\mathrm{f}}$ the fluid density, and $c_{\mathrm{p},\mathrm{f}}$ the specific heat of the fluid. The ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model lumps all radial heat transfer mechanisms into a constant effective radial thermal conductivity ${\lambda }_{\mathrm{eff},\mathrm{r}}$. The steep temperature drop at the tube walls is modeled by the introduction of an artificial wall heat transfer coefficient ${\alpha }_{\mathrm{w}}$, using Equation (2). The thermal conductivity in axial direction can be assumed to be equal to the stagnant effective thermal conductivity ${\lambda }_{\mathrm{eff},\mathrm{z}}={\lambda }_{\mathrm{eff},\mathrm{r}}^0$, or it can be neglected if the system is dominated by convective effects. The model itself has been critically discussed by many authors [31,45,46], but nevertheless widely spread due to its simplistic nature.

It was found by Yagi and Kunii [47] that the radial effective thermal conductivity can be expressed as:

$$\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}}{{\lambda }_{\mathrm{f}}}=\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}^0}{{\lambda }_{\mathrm{f}}}+\frac{1}{P{e}_{\mathrm{f},\mathrm{r}}\left(\infty \right)}PrR{e}_{\mathrm{p}}.$$

(6)

The first term on the right-hand side is the effective radial thermal conductivity of the stagnant bed. A huge number of correlations exists to determine ${\lambda }_{\mathrm{eff},\mathrm{r}}^0$, which have been reviewed by van Antwerpen et al. [48]. Based on a unit cell approach, Zehner and Schlünder [49] derived the following correlation that is widely used:

$$\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}^0}{{\lambda }_{\mathrm{f}}}=\left(1-1\sqrt{1-\epsilon }\right)+\frac{2\sqrt{1-\epsilon }}{1-{\kappa }^{-1}B}·\left(\frac{\left(1-{\kappa }^{-1}\right)B}{{\left(1-{\kappa }^{-1}B\right)}^2}ln\left(\frac{1}{{\kappa }^{-1}B}\right)-\frac{B+1}{2}-\frac{B-1}{1-{\kappa }^{-1}B}\right)$$

(7)

Here, $\kappa$ is the ratio of solid to fluid thermal conductivity and B is the deformation parameter, which is related to the void fraction by $B=1.25{\left(\left(1-\epsilon \right)/\epsilon \right)}^{10/9}$. The correlation can be further extended by incorporating secondary effects like radiative heat transfer or the effect of particle–particle contacts on the heat transfer. For a more detailed description, the interested reader is referred to the work of Tsotsas [50] and van Antwerpen et al. [48].

For the heat transfer coefficient at the wall, Yagi and Kunii [51] proposed the following correlation for $N{u}_{\mathrm{w}}={\alpha }_{\mathrm{w}}{d}_{\mathrm{p},\mathrm{v}}/{\lambda }_{\mathrm{f}}$:

$$N{u}_{\mathrm{w}}=N{u}_{\mathrm{w}}^0+\frac{1}{\left(1/N{u}_{\mathrm{w}}^{*}\right)+\left(1/N{u}_{\mathrm{m}}\right)},$$

(8)

using:

$$\begin{array}{cc}\hfill N{u}_{\mathrm{m}}& =0.054PrR{e}_{\mathrm{p}}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill N{u}_{\mathrm{w}}^{*}& =0.3P{r}^{1/3}R{e}_{\mathrm{p}}^{3/4}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill N{u}_{\mathrm{w}}^0& =\left(1.3+\frac{5}{N}\right)\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}^0}{{\lambda }_{\mathrm{f}}}.\hfill \end{array}$$

(11)

Nilles and Martin [52,53] developed the following correlation that is widely used:

$$N{u}_{\mathrm{w}}=\left(1.3+\frac{5}{N}\right)\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}^0}{{\lambda }_{\mathrm{f}}}+0.19P{r}^{1/3}R{e}_{\mathrm{p}}^{3/4}.$$

(12)

According to Dixon [31] two methods are commonly used to determine ${\lambda }_{\mathrm{eff},\mathrm{r}}$ and ${\alpha }_{\mathrm{w}}$. The first option is a parameter estimation done by conducting an optimization study based on the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, whereas the objective is to minimize the sum of squared error regarding the radial temperature profile at one or more axial positions. Alternatively, the method described by Wakau and Kaguei [54] can be used. This method is based on the approximate solution of the pseudo-homogeneous ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model and allows determining ${\lambda }_{\mathrm{eff},\mathrm{r}}$ and ${\alpha }_{\mathrm{w}}$ from the axial temperature profile in the core of the bed and the average outlet temperature. Both can easily be extracted from the particle-resolved simulation results.

The latter method was used in this study and presented in great detail by Wakao and Kaguei [54]. By neglecting the axial thermal conductivity, the analytical solution of Equation (1) is:

$$\frac{T_{\mathrm{w}}-T}{T_{\mathrm{w}}-T_0}=2\sum _{n=1}^{\infty }\left(\frac{J_0\left(2a_{n}r/D\right)exp\left(-a_n^{2}y\right)}{a_n\left(1+{\left(a_n/Bi\right)}^2\right)J_1\left(a_n\right)}\right).$$

(13)

Here, r is the radial position and $Bi$ the Biot number $Bi=\frac{{\alpha }_{\mathrm{w}}D}{2{\lambda }_{\mathrm{eff},\mathrm{r}}}$. $a_n$ is the n-th root of the following equations that include the Bessel function of the first kind and zeroth-order $J_0$ and the first kind and first-order $J_1$:

$$Bi{J}_0\left(a_n\right)=a_n{J}_1\left(a_n\right).$$

(14)

The parameter y is expressed by:

$$y=\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}z}{{\rho }_{\mathrm{f}}{u}_{\mathrm{z}}{c}_{\mathrm{p},\mathrm{f}}{\left(D/2\right)}^2},$$

(15)

whereas z is the axial position, ${\rho }_{\mathrm{f}}$ the fluid density, and $c_{\mathrm{p},\mathrm{f}}$ its specific heat. Deep in the bed, when $y\ge 0.2$, the first term in the series of Equation (13) becomes predominant, leading to:

$$\frac{T_{\mathrm{w}}-T}{T_{\mathrm{w}}-T_0}=\frac{2J_0\left(2a_{1}r/D\right)exp\left(-a_1^{2}y\right)}{a_1\left(1+{\left(a_1/Bi\right)}^2\right)J_1\left(a_1\right)}.$$

(16)

with:

$$Bi{J}_0\left(a_1\right)=a_1{J}_1\left(a_1\right).$$

(17)

In the center of the bed ($r=0$ and $T=T_{\mathrm{core}})$, Equation (16) is reduced to:

$$\frac{T_{\mathrm{w}}-T_{\mathrm{core}}\left(z\right)}{T_{\mathrm{w}}-T_0}=\frac{2exp\left(-a_1^{2}y\right)}{a_1\left(1+{\left(a_1/Bi\right)}^2\right)J_1\left(a_1\right)}.$$

(18)

If Equation (18) is logarithmized, it gives:

$$ln\left(\frac{T_{\mathrm{w}}-T_{\mathrm{core}}\left(z\right)}{T_{\mathrm{w}}-T_0}\right)=-a_1^2\left(\frac{{\lambda }_{\mathrm{eff},\mathrm{r}}}{{\rho }_{\mathrm{f}}{u}_{\mathrm{z}}{c}_{\mathrm{p},\mathrm{f}}{\left(D/2\right)}^2}\right)z+ln\left(\frac{2}{a_1\left(1+{\left(a_1/Bi\right)}^2{J}_1\left(a_1\right)\right)}\right).$$

(19)

It was shown by Wakao and Kaguei [54] that the following relationship for the average outlet temperature $T_{\mathrm{m}}$ is valid for a reasonably large axial position:

$$\frac{T_{\mathrm{w}}-T_{\mathrm{m}}}{T_{\mathrm{w}}-T_{\mathrm{core}}}=\frac{2J_1\left(a_1\right)}{a_1}.$$

(20)

From Equation (20), $a_1$ can be solved iteratively, and ${\lambda }_{\mathrm{eff},\mathrm{r}}$ can be calculated from the slope of Equation (19), subsequently. The wall heat transfer coefficient can either be determined from the intercept of Equation (19) or from Equation (17). The latter method was promoted by Wakao and Kaguei [54], since the authors argued that ${\alpha }_{\mathrm{w}}$ is very sensitive to slight changes of the intercept.

Both methods were tested during this study. A sensitivity test was conducted based on the particle-resolved CFD results for a packing of spherical particles at $R{e}_{\mathrm{p}}=100$. The sensitivity analysis of ${\alpha }_{\mathrm{w}}$ and ${\lambda }_{\mathrm{eff},\mathrm{r}}$ towards the accounted temperature range was conducted, whereas the range of ${\mathsf{\Theta }}_{\mathrm{core}}=\left(T_{\mathrm{core}}-T_0\right)/\left(T_{\mathrm{w}}-T_0\right)$ was varied as follows: $\{{\mathsf{\Theta }}_{\mathrm{core}}\in \mathbb{R}:0.05\le {\mathsf{\Theta }}_{\mathrm{core}}\le 0.4\wedge 0.6\le {\mathsf{\Theta }}_{\mathrm{core}}\le 0.95\}$. It was found that ${\lambda }_{\mathrm{eff},\mathrm{r}}$ had a relative standard deviation (RSD) of $\pm 6\%$. Then, from the intercept of Equation (19) calculated, the values of ${\alpha }_{\mathrm{w}}$ had a very low RSD of $\pm 3\%$, while the suggested method of Wakao and Kaguei increased the RSD to $\pm 15\%$, which was in contrast to the authors’ argumentation. Nevertheless, a huge discrepancy in ${\alpha }_{\mathrm{w}}$ was found: the values of the intercept-method were up to three times lower in comparison to the values determined from Equation (17). A comparison of particle-resolved CFD results against the results of the two-dimensional plug flow model in terms of axial and radial temperature profiles revealed that the temperature profiles were mispredicted if the intercept method was used, while the method promoted by Wakao and Kaguei led to reasonable results. Therefore, as Wakao and Kaguei did, we also highly recommend calculating ${\alpha }_{\mathrm{w}}$ from Equation (17) instead of the intercept of Equation (19). To evaluate ${\lambda }_{\mathrm{eff},\mathrm{r}}$ and ${\alpha }_{\mathrm{w}}$, the axial temperature profile was limited to $0.2\le \mathsf{\Theta }\le 0.8$ in this work.

2.2.2. Pseudo-Homogeneous ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$-$u_{\mathrm{z}}\left(r\right)$ Model

Instead of describing the additional thermal resistance close to the wall with a heat transfer coefficient, a radially varying effective radial thermal conductivity can be introduced. Furthermore, the radial variations of the interstitial velocity and effective axial thermal conductivity can also be considered. With this, Equation (1) is modified as follows:

$${\rho }_{\mathrm{f}}{c}_{\mathrm{p},\mathrm{f}}{u}_{\mathrm{z}}\left(r\right)\frac{\partial T}{\partial z}=\frac{1}{r}\frac{\partial }{\partial r}\left({\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)r\frac{\partial T}{\partial r}\right)+{\lambda }_{\mathrm{eff},\mathrm{z}}\left(r\right)\frac{{\partial }^{2}T}{\partial z^2}.$$

(21)

In this case, the artificial boundary condition described in Equation (2) vanishes and is replaced by the following Dirichlet boundary condition:

$$\begin{array}{ccc}\hfill T& =T_{\mathrm{w}}\hfill & \hfill \mathrm{at}r=R.\end{array}$$

(22)

As reviewed by Dixon [31], multiple models exist to determine ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$. Most often, the reactor is split into two regions to characterize the heat transfer in the near-wall and the bulk region separately. The models reported in the literature vary in their definition of the extent of each region and the description of ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)=f\left(r\right)$. Ahmed and Fahien [55] defined the wall region to be $2d_{\mathrm{p},\mathrm{v}}$ thick and used a cubic dependency for ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ in the bulk and a linear decrease in the wall region. They used the correlations of Argo and Smith [56] in combination with correlations for the radial void fraction distribution to obtain the necessary values of ${\lambda }_{\mathrm{eff},\mathrm{r}}$ in the center of the bed, at the tube wall, and at the interface of both regions. Contrary, Gunn et al. [57,58,59] used a constant value for ${\lambda }_{\mathrm{eff},\mathrm{r}}$ in the bulk region and assumed a quadratic dependency of $T\left(r\right)$ in the wall region. They defined the wall region to be $0.3d_{\mathrm{p},\mathrm{v}}$ thick. Smirnov et al. [60] defined a wall thickness that depended on bed voidage and particle specific surface area. They used a constant effective thermal conductivity in the bulk region and a linear dependency close to the wall. Winterberg et al. [61] proposed a Reynolds number-dependent thickness of the wall region. In the core region, a constant ${\lambda }_{\mathrm{eff},\mathrm{r}}$ was assumed, which decreased in the wall region, using a power-law approach that depends on the Reynolds number, Péclet number, and three more parameters. Recently, Pietschak et al. [62] reviewed several heat transfer correlations and found the correlation of Winterberg et al. [61] to be superior, especially if axial and radial variations of fluid properties were considered. Pietschak et al. [63] proposed a new correlation that accounted for the drop of ${\lambda }_{\mathrm{eff},\mathrm{r}}$ close to the wall, but without the need to introduce a discontinuity at the interface of the near wall and the bulk region. The authors correlated ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)=f\left({\rho }_{\mathrm{f}},c_{\mathrm{p},\mathrm{f}},d_{\mathrm{p},\mathrm{v}},{\epsilon }_0,{\epsilon }_{\mathrm{w}},\epsilon \left(r\right),u_0\left(r\right)\right)$ and added the cross-mixing factor and an exponent as additional parameters. The radial velocity and void fraction profiles were taken from additional correlation, but the needed data could potentially also be derived from particle-resolved simulations, as recently shown by Dixon [32].

In this work, the correlation of Winterberg et al. [61]:

$${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)={\lambda }_{\mathrm{eff},\mathrm{r}}^0+\frac{u_{0,\mathrm{c}}{d}_{\mathrm{p},\mathrm{v}}{\rho }_{\mathrm{f}}{c}_{\mathrm{p},\mathrm{f}}}{K_{\mathrm{w}}}·f\left(R-r\right),$$

(23)

using:

$$f\left(R-r\right)=\left\{\begin{array}{cc}{\left(\frac{R-r}{k_{\mathrm{f},\mathrm{w}}{d}_{\mathrm{p},\mathrm{v}}}\right)}^{n_{\mathrm{f},\mathrm{s}}}\hfill & \mathrm{if}0<R-r<k_{\mathrm{f},\mathrm{w}}{d}_{\mathrm{p},\mathrm{v}}\hfill \\ 1\hfill & \mathrm{if}{k}_{\mathrm{f},\mathrm{w}}{d}_{\mathrm{p},\mathrm{v}}<R-r<R,\hfill \end{array}\right.$$

(24)

was used as the basis to determine ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$. In Equation (23), $K_{\mathrm{w}}$ is the cross-mixing factor that describes the relationship between effective thermal conductivity, particle shape, and flow velocity deep in the bed. The cut-off parameter $k_{\mathrm{f},\mathrm{w}}$ in Equation (24) sets the dimensionless wall distance, after which the constant thermal conductivity, which was assumed in the core region, drops towards the wall. The exponent $n_{\mathrm{f},\mathrm{s}}$ describes the curvature of the damping function close to the wall.

To determine the parameters above, the circumferentially averaged radial temperature profiles in the interval $z=\left[0.1:0.1:1.1\right]$ were extracted from the particle-resolved CFD simulations for the simulation with $R{e}_{\mathrm{p}}\ge 500$. For the lowest investigated Reynolds number of $R{e}_{\mathrm{p}}=100$, the radial temperature gradients flattened out quickly. Therefore, in this case, only the temperature profiles in the range of $z=\left[0.1:0.1:0.3\right]$ were considered. Based on the model described by Equation (21), a parameter optimization study was conducted, using the Nelder–Mead algorithm, while the objective was to minimize the sum of least squares of the difference of the radial temperature profiles between the simplified model and the particle-resolved results. In total, a number of 1100 ($R{e}_{\mathrm{p}}\ge 500$) and 300 ($R{e}_{\mathrm{p}}=100$) data points were available for the optimization task for each operating condition.

3. Heat Transfer Validation

Experimental validation data for axial or radial temperature profiles are scarce and hard to find. Nevertheless, Wehinger et al. [3] and Dong et al. [6] were able to prove the accuracy of the particle-resolved CFD approach, especially in combination with the local “caps” meshing strategy, in terms of axial and radial temperature profiles.

Based on experimental data that were provided by Clariant International Ltd., a validation study was conducted in this work to confirm the reliability of particle-resolved CFD also under industrially relevant conditions (T > 1000 °C). The experimental setup consisted of a hot box, fired with an electrical furnace, and a single reformer tube with an inner diameter of 0.1016 m and a bed height of 1 m. With thermocouples, placed at the outside of the reformer tube, the axial profile of the outer wall temperature was measured. The temperature in the center of one of the packed reformer tubes was measured with a 0.25″ standard 316 SS axial thermowell until an axial distance of approximately 0.5 m. The thermocouples used were of type K, with an accuracy of ±1.5 °C or ±0.4%, whichever was greater. In the scope of this study, two different particle shapes, a tablet-like cylinder with six holes (33 × 18 mm) and an almost equilateral cylinder with ten holes (19 × 16 mm), were investigated.

The experimental setup was replicated numerically, whereas special emphasis was given to meet the particle count that was determined in the experiments. To achieve this, particle static friction coefficients were calibrated, as described by Jurtz et al. [4,42]. Since the tube thickness was relatively big, not only the fluid and the particles, but also the reformer tube itself were spatially discretized. While a fully conformal contact interface was used for the particle-fluid interface, for the sake of reduced cell count, the tube-fluid interface was performed as a non-conformal mapped contact interface. An overview of the investigated setups, including snippets of the resulting meshes, can be seen in Figure 3. Preliminary studies showed that the thermowell not only affected the flow field significantly, as recently discussed by Dixon and Wu [64]. It was found that the heat conduction through the thermowell could not be neglected, since it significantly affected the temperature distribution in the vicinity of the measuring device. To consider for heat conduction in the thermowell, a three-dimensional shell model was used that solved for the lateral conductive heat transport and modeled the heat transport in the face normal direction via the assumption of a constant temperature gradient. The radiative heat transport was considered using a surface-to-surface radiation model as done by Wehinger at al. [7] recently.

The experimentally measured temperature distribution on the outer side of the reformer tubes was applied as a spatially varying fixed temperature boundary condition at the outer tube surface. The inlet temperature, according to experimental data, was set to 560 °C and the operating pressure to 1.5 barg. Nitrogen was used as the working fluid, whereas the ideal gas equation of state was used. The fluid viscosity and thermal conductivity were calculated using the Chapman–Enskog model. Inlet flow rates were varied between 15 and 50 Nm³/h. Particles’ thermal conductivity was set to 0.25 W/(m K), whereas for the tube and thermowell, the following function, derived from the spec sheet, was used: ${\lambda }_{\mathrm{s}}\left[W/\left(mK\right)\right]=8.195·exp\left(1.188·{10}^{-3}·T\right)$. The emissivity was set to 0.75 for the particles surface and to 0.6 for the inner reformer tube and the thermowell.

The simulation results are given in Figure 4 in terms of the axial profiles of the dimensionless temperature $\mathsf{\Theta }=\left(T-T_0\right)/\left(T_{\mathrm{ref}}-T_0\right)$, whereas $T_{\mathrm{ref}}$ is the furnace temperature. The numerical data are presented as a scattered cloud of small symbols to also visualize the temperature variation in the circumferential direction. It can be seen that due to the conductive heat transport within the solid of the thermowell, temperature variations in circumferential direction were low. Without considering this heat transfer mechanism, temperature differences of over 50 K were found (see Section S3), which indicated that the measuring device not only affected the fluid dynamics, but also the measured temperature field significantly. This strengthened the argument that the use of high-fidelity numerical methods can improve the accuracy of determining effective heat transfer parameters significantly. Deep in the bed, an excellent agreement could be found between the predicted and measured temperatures for the six-hole tablets. Only for $z/h\approx 0.1$, some deviations were found. However, if one considers the obvious impact of the heterogeneous bed morphology on the axial temperature profile, the accuracy was still acceptable. For the 10-hole cylinders, the experimental temperature profile was hit almost perfectly for $z/h\le 0.35$. For a flow rate of 15 N m³/h, also deep in the bed, a good agreement with the experimental data was found. However, at higher flow rates, for $z/h\approx 0.5$, the simulations results were far off. The reason for this could not be identified, but the fact that the experimental data showed an increase of the axial temperature gradient at higher bed depths for high flow rates was suspicious and may indicate that the temperature sensors were damaged under the harsh operating conditions. Similar problems have been noted by other authors [31] and illustrate the challenges associated with experimental temperature measurements in fixed beds. A possible re-ordering of particles at the tip of the thermowell during operation might also be a possible reason for the deviations observed.

4. Results and Discussion

The heat transport in fixed-bed reactors is strongly coupled to fluid dynamics effects that are induced by the heterogeneous bed morphology. Therefore, in the first part, the bed morphology and fluid dynamics of all generated packings were investigated. Afterwards, the different configurations were compared with regard to their thermal performance. The global heat transfer coefficient $U={\dot{Q}}_{\mathrm{w}}/\left(A_{\mathrm{w}}{\mathsf{\Delta }}_{\log }T\right)$, using ${\mathsf{\Delta }}_{\log }T=\left(T_{\mathrm{out}}-T_{\mathrm{in}}\right)/log\left(\left(T_{\mathrm{w}}-T_{\mathrm{in}}\right)/\left(T_{\mathrm{w}}-T_{\mathrm{out}}\right)\right)$, was used to compare the different designs, whereby U was evaluated for an axial threshold of the reactor that fulfilled the criterion $0.0\le {\mathsf{\Theta }}_{\mathrm{core}}\le 0.8$. In the last part, the simulation results were used to determine effective thermal transport parameters that are commonly needed for the pseudo-homogenous two-dimensional plug-flow model. The results were compared against particle-resolved CFD results to understand the reliability of simplified models. A summary of the most important results and simulation parameters is given in

Table 2. Summary of important results and simulation parameters.

Shape		General Properties					${\mathit{\lambda }}_{\mathbf{eff}\mathbf{,}\mathbf{r}}\mathbf{\text{-}}{\mathit{\alpha }}_{\mathbf{w}}\mathbf{\text{-}}\mathbf{Model}\mathbf{Parameters}$		${\mathit{\lambda }}_{\mathbf{eff}\mathbf{,}\mathbf{r}}\left(\mathit{r}\right)\mathbf{\text{-}}\mathbf{Model}\mathbf{Parameters}$
		$\mathit{R}\mathit{e}\mathbf{\left(}\mathbf{\text{-}}\mathbf{\right)}$	$\mathbf{\epsilon }\mathbf{\left(}\mathbf{\text{-}}\mathbf{\right)}$	$\mathsf{\Delta }\mathit{p}\mathbf{/}\mathsf{\Delta }\mathit{z}\mathbf{(}\mathbf{Pa}\mathbf{/}\mathbf{m})$	${\mathit{\lambda }}_{\mathbf{eff}\mathbf{,}\mathbf{r}}^{\mathbf{0}}\mathbf{(}\mathbf{W}\mathbf{/}\mathbf{(}\mathbf{m}\mathbf{K}\mathbf{\left)}\mathbf{\right)}$	U (W/(m2 K))	αw (W/(m2 K))	${\mathit{\lambda }}_{\mathbf{eff}\mathbf{,}\mathbf{r}}\mathbf{(}\mathbf{W}\mathbf{/}\mathbf{(}\mathbf{m}\mathbf{K}\mathbf{\left)}\mathbf{\right)}$	${\mathit{K}}_{\mathbf{n}}\mathbf{\left(}\mathbf{\text{-}}\mathbf{\right)}$	${\mathit{n}}_{\mathbf{fw}}\mathbf{\left(}\mathbf{\text{-}}\mathbf{\right)}$	${\mathit{k}}_{\mathbf{fw}}\mathbf{\left(}\mathbf{\text{-}}\mathbf{\right)}$	MSE (K2)
Spheres	loose	100	0.473	50.0	0.07943	30.86	108.44	0.4048	7.487	1.017	0.2808	1.08
		500	0.473	651.2	0.07943	44.73	172.74	1.4251	10.210	1.242	0.3210	1.40
		1000	0.473	2206.5	0.07943	59.08	236.83	2.3836	10.741	1.246	0.3431	2.41
		2000	0.473	8429.9	0.07943	95.09	369.63	4.9837	11.448	1.310	0.3031	3.26
	dense	100	0.426	65.0	0.08759	20.55	255.61	0.2726	12.051	0.566	0.2711	12.76
		500	0.426	837.9	0.08759	26.83	192.19	1.0301	18.974	0.955	0.3290	5.78
		1000	0.426	2868.8	0.08759	40.58	276.36	1.8464	20.127	1.004	0.3290	6.22
		2000	0.426	11,258.8	0.08759	67.95	466.10	3.7963	20.111	1.012	0.3325	7.65
Cylinders	loose	100	0.463	76.3	0.10052	49.82	88.92	0.5529	5.308	1.253	0.3376	0.84
		500	0.463	1021.5	0.10052	67.83	161.63	2.0920	8.956	1.995	0.1574	3.40
		1000	0.463	3539.1	0.10052	94.47	241.29	3.8584	8.855	1.981	0.1624	2.36
		2000	0.463	14,417.8	0.10052	151.31	408.46	7.9587	8.931	2.027	0.1430	2.66
	dense	100	0.371	158.1	0.12021	46.08	180.10	0.4756	6.042	1.536	0.2753	5.00
		500	0.371	2059.1	0.12021	58.19	279.40	1.9381	9.652	1.448	0.2477	3.48
		1000	0.371	7254.7	0.12021	88.99	420.27	3.5749	10.697	1.532	0.2087	3.89
		2000	0.371	33,613.3	0.12021	151.82	756.64	7.8851	10.721	1.478	0.2019	4.29
Rings	loose	100	0.755	37.6	0.05826	36.40	45.92	0.2522	5.692	1.030	0.5181	2.94
		500	0.755	504.2	0.05826	64.91	86.39	1.4277	5.794	1.435	0.3768	0.97
		1000	0.755	1743.2	0.05826	88.33	126.61	2.9977	5.815	1.605	0.3218	1.24
		2000	0.755	6583.6	0.05826	124.33	198.75	5.9213	5.816	1.686	0.2936	1.65
	dense	100	0.710	65.9	0.06587	32.18	49.04	0.3024	4.125	0.653	1.8197	4.46
		500	0.710	891.0	0.06587	69.53	99.14	1.3419	4.177	1.022	0.5710	2.59
		1000	0.710	3136.8	0.06587	95.12	157.84	2.3876	6.289	1.285	0.4061	2.01
		2000	0.710	12,675.6	0.06587	146.97	257.85	5.4581	6.324	1.492	0.3022	1.91
4-hole cylinders	loose	100	0.634	55.9	0.07945	33.10	68.00	0.2902	7.492	1.199	0.3727	1.93
		500	0.634	671.3	0.07945	60.05	117.27	1.3121	7.243	1.403	0.3596	1.15
		1000	0.634	2271.6	0.07945	87.79	167.40	2.8781	7.012	1.518	0.3252	1.60
		2000	0.634	8622.2	0.07945	139.90	275.54	5.9914	6.963	1.634	0.2744	2.18
	dense	100	0.571	105.7	0.09141	29.57	113.27	0.2724	5.641	0.604	1.4786	3.93
		500	0.571	1257.3	0.09141	54.84	143.32	1.2905	6.981	0.912	0.7643	2.71
		1000	0.571	4327.0	0.09141	87.87	182.08	3.1602	7.514	1.107	0.5110	3.08
		2000	0.571	17,790.1	0.09141	146.86	351.64	6.0885	6.748	1.111	0.5690	3.50
Wall structure	loose	100	0.496	41.7	0.07570	38.08	101.09	0.4923	7.065	0.424	0.5479	3.16
		500	0.496	560.4	0.07570	61.63	176.84	1.9282	8.502	0.839	0.4355	4.23
		1000	0.496	1934.4	0.07570	85.03	317.63	2.9432	9.177	1.073	0.3126	4.71
		2000	0.496	7640.2	0.07570	131.56	544.53	5.6290	9.806	1.308	0.2195	4.87
	dense	100	0.439	67.6	0.08525	38.32	252.46	0.4644	7.057	1.121	0.1120	1.57
		500	0.439	888.3	0.08525	61.54	308.36	1.5457	9.830	1.098	0.2072	1.56
		1000	0.439	3078.3	0.08525	87.89	424.75	2.9964	10.592	1.210	0.1952	2.12
		2000	0.439	12,539.7	0.08525	149.93	713.48	6.5500	11.011	1.416	0.1525	2.57

. A schematic drawing of the setup is presented in Figure 5.

4.1. Bed Morphology and Fluid Dynamics

Already, the first visual impression of the generated packings that was given in Figure 1 showed the strong impact that the packing mode had on the particle arrangement. This can best be seen for spherical and cylindrical particles. While for the loose packing configuration, although the confining walls exerted an ordering effect on the particles, to some extent, random arrangement of the particles close to the wall can be seen, the compacted beds were characterized by a high degree of order. Especially the spherical particles tended to build band-like structures at the wall, whereas cylindrical particles built stacked structures and were mostly oriented parallel or perpendicular to the wall. From a fluid dynamics and reaction engineering point of view, the most important effect was the significant reduction in bed voidage that was caused by bed densification. By this, the pressure drop, local flow phenomena, hydraulic residence time, and the active catalytic surface area per reactor volume were significantly affected. The evaluated bed voidage, listed in Table 2, shows that for spherical particles, the bed voidage was reduced by 10%. An extreme reduction of 20% was found for cylindrical particles. For particles with inner voids, like rings and four-hole cylinders, the effect was less pronounced, giving a drop of 10% and 6%, respectively. However, this reduced impact was only a result of the overall higher bed voidage for these particles. For the configuration of spherical particles in the reactor with macroscopic random wall structures, the densification-induced reduction of bed voidage was 11%, which was similar to the reactor with plain walls.

The axial and radial void fraction profiles were good resources to understand the packing morphology of the different designs. Strong and regular oscillations are indicators of ordered particle arrangements and a loss of randomness in the system, whereas low non-regular fluctuations in bed voidage point towards an increasing randomness of the particle arrangement. For an ideally random packing arrangement, the void fraction profile should end in a constant value. Distinct peaks in the void fraction profile are indicators of additional voids that are a result of a non-appropriate filling strategy, which leads to jamming of particles. The axial void fraction profiles of all investigated packings are given in Figure 6. Since the bed rested on a bottom plate, the lowest layers of particles experienced a certain ordering effect, which was induced by the adjacent wall. For spherical particles, only a point contact was possible between the particles and the bottom plate, leading to a value of $\epsilon =1$ at $z/d_{\mathrm{p}}=0$. Particles of the cylindrical shape type may have a point, line, or face contact with the wall. If face contacts are present, it is possible that $\epsilon <1$ at $z/d_{\mathrm{p}}=0$. However, for most of the investigated packing, it can be seen that the ordering effect of the bottom plate led to regular oscillations in the void fraction that flattened out after a distance of 3–d_p and ended up in random oscillations of lower magnitude, indicating a stochastic axial distribution of the particles. The only exceptions were the compacted packing of spherical particles and the loose packing of spheres in the reactor with macroscopic wall structures. For the dense packing of spheres, regular oscillations were observed between $0\le z/d_{\mathrm{p}}\le 22$. This indicated a pronounced layer formation in the bottom part of the reactor. In the remaining part of the reactor as well, regular oscillations were observed, albeit to a lesser extent. In the wall structured reactor, high fluctuations were observed that suddenly appeared and flattened out. A probable reason for this was jamming of particles during the filling process that led to additional voids. This hypothesis was strengthened by the fact that this effect vanished for the densified packing.

Of fundamental interest for the understanding of the fluid dynamics are the radial void fraction profiles and the radial profiles of the circumferentially averaged axial velocity, given in Figure 7. Here, the axial velocity was normalized to the local interstitial velocity $u_0/\epsilon \left(r\right)$. With the exception of the structured wall reactor, for all particle shapes, directly at the wall, a void fraction of $\epsilon =1$ was found due to the presence of point and/or line contacts, only. For spherical particles, a first minimum in the void fraction was reached after the distance of one particle radius away from the wall, indicating that the majority of spheres were in direct contact with the wall, forming a closed particle layer. Furthermore, local minima and maxima occurred at positions corresponding to multiples of the particle radius, whereas the oscillations slightly decreased. The global minimum of the void fraction was located in the center of the bed, indicating that an almost stacked arrangement of spheres was present. This was the result of odd tube-to-particle diameter ratios [40,41]. A strong correlation could be found between the void fraction and the velocity profile. Close to the wall, the velocity reached its maximum, known as the wall channeling effect. The position of further minima and maxima corresponded directly to the position where high/low void fractions were found. While the minima of axial velocity did not change with varying $R{e}_{\mathrm{p}}$, the maxima increased slightly in the center of the bed if $R{e}_{\mathrm{p}}$ was lowered. This effect could be attributed to the gas expansion due to heating and to the decreasing wall effect if $R{e}_{\mathrm{p}}$ was lowered. The center of the bed was almost completely blocked for the flow. The above findings were also valid for the densified packing of spheres; however, the effects were even more pronounced, resulting in a complete blockage of flow paths at $r^{*}=\left(R-r\right)/d_{\mathrm{p}}=\left[0.5,1.5,2.5\right]$, and strong channeling was observed at $r^{*}=\left[0.1,1.0,2.0\right]$, whereas for $R{e}_{\mathrm{p}}\le 500$, the strongest channeling was not found at the wall, but at $r^{*}=2.0$, which is very uncommon.

For cylindrical particles, the trend was similar as for spheres; however, the minima/ maxima in the void fraction and velocity were slightly shifted towards the bed center, which indicated that some particles were diagonally aligned. For the densified packing, the minima/maxima were found at multiples of the particle radius, which was a result of the particles’ preferred parallel/orthogonal alignment. In contrast to the packings of spheres, where the wall channeling was almost independent of $R{e}_{\mathrm{p}}$, for cylindrical particles, the wall channeling effect increased significantly if $R{e}_{\mathrm{p}}$ was raised. This effect became very dominant for the compacted packing. The void fraction profiles for Raschig rings and four-hole cylinders looked pretty complex; nevertheless, especially for $r^{*}<1$, the inner voids of the particles were clearly reflected by corresponding additional maxima in the void fraction. However, no maxima in velocity could be found at void fraction maxima that corresponded to inner voids. This indicated that the flow through the inner particle voids was partially blocked, which might be because of an orthogonal particle alignment. Overall, the void fraction and velocity oscillations were less pronounced for those particle shapes, but heterogeneities increased if the beds were compacted. Similar to cylindrical particles, the wall channeling effect increased with $R{e}_{\mathrm{p}}$ and became more pronounced for densified packings.

The use of macroscopic random wall structures for packings of spherical particles changed the void fraction and velocity profiles significantly. Due to the presence of the wall structure, the void fraction at the wall fell to a value of $\epsilon \approx 0.56$. As a result, the wall channeling effect was hindered, and fluctuations in the void fraction and velocity were qualitatively more comparable to the ones of Raschig rings than spheres. The densification of the bed led to slightly more pronounced minima and maxima; however, this effect was not as distinct as for spherical particles in a smooth walled reactor.

4.2. Heat Transfer Characteristic

A fair comparison of the thermal performance of different reactor concepts always depends on the process boundary conditions that are set. Figure 8 shows the global heat transfer coefficient U as a function of different parameters. Re-fitting of an existing unit that is integrated in a complex production process can lead to the necessity of keeping the throughput constant, which is equivalent to keeping $R{e}_{\mathrm{p}}$ invariant. In this case, especially at low $R{e}_{\mathrm{p}}$, cylindrical particles showed the most beneficial heat transfer characteristic, followed by the wall-structured reactor, Raschig rings, and four-hole cylinders. Spherical particles performed worst over the complete range of investigated $R{e}_{\mathrm{p}}$. At high $R{e}_{\mathrm{p}}$, cylindrical particles still performed best; however, rings, four-hole cylinders, and the reactor with wall structures were close. For spheres, cylinders, and four-hole cylinders, U decreased if the packings were compacted. This is of special interest, since in industrial applications, most often, densified packings are used to ensure the same pressure drop in the different tubes of the tube bundle reactor. Interestingly the effect was less pronounced for Raschig rings and the wall-structured reactor. At high $R{e}_{\mathrm{p}}$, even a slight increase in thermal performance could be seen for those reactor types. In general, the performance gain induced by macroscopic random wall structures was significant.

Another valid process boundary condition can be the necessity of keeping the hydraulic residence time invariant. In this case, $R{e}_{\mathrm{p}}/\epsilon$ needs to be kept constant. Under this constraint, Raschig rings and four-hole cylinders performed best for moderate to high $R{e}_{\mathrm{p}}$, followed by cylinders and the wall-structured reactor. For the lowest investigated $R{e}_{\mathrm{p}}$, again, cylindrical particles seemed to perform slightly better than rings.

If a new plant is built and process-driven constraints are low, the most energy efficient particle shape might be an appropriate choice. In this case, the specific pressure drop $\mathsf{\Delta }p/\mathsf{\Delta }z$ can be one parameter that should be kept constant when comparing different designs. In this case, Raschig rings, the wall-structured reactor, and cylinders performed best. The comparison of the designs from this energetic point of view showed that bed densification led to a less energy-efficient thermal performance, whereas this effect was less pronounced for Raschig rings and the reactor with macroscopic wall structures.

4.3. Effective Thermal Transport Properties

As discussed, particle-resolved CFD is a valuable tool to support process intensification on the meso-scale level, e.g., by finding optimized particle shapes [4,10,11,12] or new reactor concepts, e.g., by applying macroscopic wall structures [20,21] or using internals [19]. However, for process intensification on a macroscopic scale, e.g., by running plants under dynamic operation conditions, developing process integration strategies, or doing plant optimization, different numerical tools are necessary. Process simulation platforms often use pseudo-homogeneous two-dimensional plug flow models. Depending on the class of model used, certain effective transport parameters are needed, which are often not known. In this section, methods are presented for how those parameters can be extracted from particle-resolved CFD results.

4.3.1. ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ Model

Although its limitations are well known [31], the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model is still widely spread, due to its efficiency and simple implementation. Here, the radial heat transport was characterized by the wall heat transfer coefficient ${\alpha }_{\mathrm{w}}$ and the effective radial thermal conductivity ${\lambda }_{\mathrm{eff},\mathrm{r}}$, which was assumed to be uniform everywhere in the reactor. By extracting the axial core temperature profile and average inlet/outlet temperatures from the CFD simulations, both parameters were determined by using Equations (17), (19), and (20). The results are summarized in Table 2. The parameters were then used to calculate the temperature fields by using the pseudo-homogeneous model described by Equation (1) in conjunction with the boundary condition in Equation (2).

A one-to-one comparison of all investigated cases in terms of radial temperature profiles at different axial positions is provided in Supplementary Material, Section S4. A condensed visualization of the results is given in Figure 9. Here, the deviations of the circumferentially averaged temperature fields, predicted by the pseudo-homogeneous model, are given in relation to the particle-resolved CFD results. Deep red and deep blue colors indicate that the deviation was above or below 10 K. This critical cut-off temperature was chosen, motivated by the rule of van’t Hoff, saying that the speed of a chemical reaction doubles to triples itself when the temperature is raised by 10 K [65]. The characteristic temperature drop at the wall that was a result of ${\alpha }_{\mathrm{w}}$ can only hardly be seen in Figure 9. The reader is referred to the radial temperature profiles given in Section S4. It can be seen that the temperatures close to the wall (r^* ≤ 0.2–0.4) were systematically underpredicted by the simplified model. This drawback is well known and deeply discussed by many authors [31]. Furthermore, the model was not able to capture morphological and fluid dynamic heterogeneities, which led to step-like temperature profiles, as can be seen best for the radial temperature profiles of the dense spherical packing. Recently, this was also found by Moghaddam et al. [33], who introduced heterogeneities by increasing the solid thermal conductivity. Besides those systematic errors, the deviation in relation to the CFD results was relatively low for the majority of cases. For all investigated designs, the deviation was less than 5 K for the rear part of the reactor ($z/d_{\mathrm{p}}\ge 40$). The threshold of 10 K was mostly exceeded in the entry zone ($z/d_{\mathrm{p}}\le 20$). Overall, there seemed to be a trend that deviations increased if $R{e}_{\mathrm{p}}$ was raised. The method seemed to work equally well for loose and densified packings with a slight trend towards less deviations for dense beds. Considering the numerical effort that the simplified model needed in comparison to the particle-resolved CFD simulation, which was ≈10 s compared to ≈24 h, the accuracy was still remarkable.

4.3.2. ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$-Model

The obvious drawback of the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, that the additional near wall thermal resistance is only captured with an artificial temperature drop directly at the reactor wall, can be circumvented by using a radially varying effective radial thermal conductivity. In this work, the correlation of Winterberg [61] (see Equation (23)) was used as the basis to determine the effective radial thermal conductivity. The three necessary parameters of the Winterberg correlation were determined by conducting a parameter optimization study. The basis of this study was the transport equation described by Equation (21). A summary of the optimized model parameters, including the mean squared error $MSE$, is given in Table 2.

The comparison of the radial temperature profiles for different axial positions can be found in Supplementary Material, Section S5. The spatially resolved deviations between the simplified model and the CFD results are given in Figure 10. It is obvious that in comparison to the results of the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, the accuracy was significantly improved. The temperature close to the reactor wall was predicted with a high degree of accuracy by the model. Only sporadically, the temperatures were overestimated by more than 10 K in the vicinity of the wall, whereby the location was mostly limited to the entry zone ($z/d_{\mathrm{p}}\le 20$). A direct comparison of the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model and the ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ model in relation to the CFD results is given in Figure 11 for the loose packing of Raschig rings at $R{e}_{\mathrm{p}}=1000$, showing the superior accuracy of the ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ model, especially close to the wall. Deep in the bed, deviations outside of the 10 K threshold were mostly found for packings that were characterized by a higher degree of morphological heterogeneity, like the packings of cylindrical particles, the dense bed of spheres, and the loose packing of spheres in the reactor with a random wall structure. While the former configurations were characterized by strong variations in the radial void fraction distribution, the latter showed big fluctuations in the axial void fraction profile. Bigger deviations were mostly limited to the entry zone, indicating that thermal entrance effects, which were not resolved by an axially invariant ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$, might be the reason for this. In contrast to the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, the deviations did not seem to increase if $R{e}_{\mathrm{p}}$ was raised. Since the Winterberg correlation did not explicitly consider local variations in the void fraction or axial velocity, it was, similar to the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, not able to capture the step-like effects of the temperature profiles.

5. Conclusions

In this work, it was shown in which way the particle-resolved simulation of fixed-bed reactors can play a central role in the process the intensification of this reactor type. After a brief validation study, showing that also under harsh industrial conditions, particle-resolved CFD was able to predict the temperature field accurately, the heat transfer characteristics of different particle designs were investigated. The studied designs differed in the used particle shape and the bed density. The results showed that heterogeneities in radial void fraction distribution and axial velocity increased, if packings were compacted. As a result, the overall heat transfer coefficient U decreased for most particle shapes. Although the wall channeling effect was most pronounced for the fixed-beds of cylindrical particles, it was found that this particle shape was among the most efficient, with respect to U. Furthermore, a novel reactor tube design that used random macroscopic wall structures was investigated. For packings of spherical particles, it was found that macroscopic random wall structures can significantly decrease morphological heterogeneities, leading to a significantly better heat transfer characteristic. Taking into account various process-related boundary conditions, cylindrical particles, Raschig rings, and wall-structured reactors were identified as the most promising concepts to intensify the radial heat transport.

Methods were presented to determine effective thermal transport parameters, which are needed for simplified pseudo-homogeneous models, from the particle-resolved CFD results. Depending on the degree of morphologically induced heterogeneities, an excellent to fair agreement was found for the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model in comparison to the CFD results, whereas deviations became bigger if the morphology was more heterogeneous. The known problem of underestimated temperatures close to the reactor tube, as one of the biggest drawbacks of this model, was confirmed. To circumvent this problem, parameter optimization studies, based on the Winterberg correlation, were performed to predict the radially varying effective radial thermal conductivity, which was needed for the ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ model. A very good agreement regarding the radial temperature profiles was found between the ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ model and the particle-resolved CFD results. In comparison to the ${\lambda }_{\mathrm{eff},\mathrm{r}}$-${\alpha }_{\mathrm{w}}$ model, the ${\lambda }_{\mathrm{eff},\mathrm{r}}\left(r\right)$ model showed its superior accuracy close to the reactor wall. Nevertheless, it was found that both pseudo-homogeneous models became less accurate if step-like temperature profiles, which were a either a result of morphological heterogeneities or a high particle solid thermal conductivity, were present.

In terms of process intensification, this work showed that particle-resolved CFD can either directly be used to study improvements on the meso-scale through:

• studying the impact of particle shape, internals, or reactor tube design on the performance;
• investigating the effect of operating conditions and physical properties;
• testing of novel reactor tube concepts, e.g., reactors with random macroscopic wall structures or heat fins;
• identifying local phenomena as hot/cold spot formation or catalyst poisoning;

or as a reliable source for parameters, and correlations of those, that are needed for process simulation. This allows a more reliable analysis of process intensification by studying:

• dynamic operating conditions;
• process integration concepts;
• conducting process design optimization.