Influence of Macroscopic Wall Structures on the Fluid Flow and Heat Transfer in Fixed Bed Reactors with Small Tube to Particle Diameter Ratio

Abstract

Fixed bed reactors are widely used in the chemical, nuclear and process industry. Due to the solid particle arrangement and its resulting non-homogeneous radial void fraction distribution, the heat transfer of this reactor type is inhibited, especially for fixed bed reactors with a small tube to particle diameter ratio. This work shows that, based on three-dimensional particle-resolved discrete element method (DEM) computational fluid dynamics (CFD) simulations, it is possible to reduce the maldistribution of mono-dispersed spherical particles near the reactor wall by the use of macroscopic wall structures. As a result, the lateral convection is significantly increased leading to a better radial heat transfer. This is investigated for different macroscopic wall structures, different air flow rates (Reynolds number Re = 16 ...16,000) and a variation of tube to particle diameter ratios (2.8, 4.8, 6.8, 8.8). An increase of the radial velocity of up to 40%, a reduction of the thermal entry length of 66% and an overall heat transfer increase of up to 120% are found.

1. Introduction

Since the early 1950s, numerous authors have investigated fixed bed reactors regarding the packing structure [1], fluid flow [2] and heat and/or mass transfer [3]. Several mathematical models have been developed to predict the pressure drop and to describe heat and mass transfer phenomena. These models can be divided into three fundamentally different categories:

• Pseudo-homogeneous models, in which particles are not explicitly resolved. Their impact on pressure drop and heat/mass transfer is considered by using effective transport parameters that incorporate all non-resolved phenomena. Recently, several heat transfer correlations for pseudo-homogeneous reactor modeling were evaluated [4] and found, that most accurate results can be achieved if the effective radial thermal conductivity is calculated using the correlation of Winterberg and Tsotsas [5].
• Heterogeneous models, e.g., the model developed by Schluender and Tsotsas, where different sets of equations are solved for the fluid and the solid phase. Nevertheless, the pellets are not geometrically resolved, and the coupling between the fluid and the solid phase relies on empirical correlations for heat and mass transfer parameters.
• Particle-resolved computational fluid dynamics (CFD), where each pellet and the fluid, including the interstices between the particles are fully spatially resolved, and Navier–Stokes equations are solved. A detailed description of this modeling approach is presented in the comprehensive review articles by Dixon [6], Jurtz et al. [7] and Dixon and Partopour [8].

It was shown by Roblee et al. [1] that for fixed beds with a large tube to particle diameter ratio, the assumption of a constant void fraction and a plug flow profile suits quite well. On the other hand, those assumptions are not valid in multi-tubular packed bed reactors, which are widely used in the chemical and process industry for highly exothermic or endothermic catalytic reactions. Multi-tubular packed bed reactors consist of several thousands of tubes with a small tube to particle diameter ratio. In those reactor tubes a strong coupling exists between the local packing structure, velocity and temperature field. The higher voidage near the reactor wall leads to an increase of axial velocity in this region. This effect is well known as the bypass effect and inhibits the convective heat transport from the bed to the reactor wall or vice versa. A good heat transfer is important due to safety and quality requirements, and an improvement of the radial heat transfer is aspired.

Although the complex characterization of fluid flow and heat transfer in fixed beds with small tube to particle diameter ratio is not yet finally concluded it is well known that the bypass effect has a non-negligible impact on the heat transport near the reactor wall. According to Yagi and Kunii [9], the increasing heat transport resistance near the wall region is caused by the following effects:

• decreasing impact of the solid conductivity due to the higher voidage near the wall;
• development of a laminar boundary layer;
• decreasing lateral mixing due to the bypass effect.

A previous study by Zobel et al. [10] has shown that the use of a macroscopic wall structure leads to a more homogeneous radial void fraction distribution, and as a consequence, to a higher averaged radial velocity. It was concluded that this higher radial velocity leads to a better lateral mixing, but no evidence is given and further investigations are proposed to quantify the effect and to find an optimal wall structure. To the best of our knowledge, no further investigations in that direction were conducted or published. Particle-resolved 3D simulations are a cost-effective way to simulate and optimize or at least investigate different designs due to its relatively easy setup without any physical construction and its detailed and easy access to flow field information.

Particle-resolved simulation of fixed bed reactors is coupled with two major challenges: Firstly, a CAD description of a representative bed of particles is required. For this purpose, several authors used the statistical approach of the Monte Carlo Method [11,12,13] while others took advantage of the discrete element method (DEM), which is physically motivated and replicates the filling process [10,14,15]. Secondly a suitable mesh must be generated which avoids bad quality cells in the vicinity of contact points between particles and between the particles and the confining wall. Different approaches have been developed during the last years, which can be classified into global and local methods. It is known that the global methods like shrinking or inflating of the particles have a significant impact on the bed voidage and therefore on the calculated pressure drop [14]. This impact can be reduced by a meshing strategy which modifies the packing only locally like bridging two adjacent particles with a small cylinder [16] or by flattening the particles locally at their contact point [15].

In this work the effect of different macroscopic wall structures on the radial and axial void fraction and flow distribution was investigated. The impact of the Reynolds number Re and the tube to particle diameter ratio $D/d_p$ is quantified. Furthermore, the question is addressed whether a higher radial velocity has a positive influence on the radial convective heat transport and therefore on the thermal characteristic of the reactor. The results clearly show that with appropriate modifications of the walls, a significantly higher heat transfer can be achieved, opening the possibility for new operating conditions.

2. Numerical Methods

In this work, DEM is used for the generation of the fixed bed while the local flattening meshing strategy [15] is adopted for the mesh generation with some slight modifications. These modifications affect only the workflow, so in contrast to using three different software packages as in the original publication, only Simcenter STAR-CCM+ is now used for the DEM, CAD preparation, meshing and flow calculation in a fully automated way. The calculation of the fluid flow and the temperature field is based on the finite volume method.

2.1. DEM

To generate a random packing, mono-dispersed spherical particles are initialized randomly in the upper part of the tube and fall because of gravity to the bottom of the tube. For each of these particles a force balance is formulated and solved, which takes gravity, the interaction between the particles and the interaction between the particles and the tube walls into account:

$$m_i\frac{d{v}_i}{dt}=F_g+\sum _{np}{F}_c+\sum _{nb}{F}_c$$

(1)

and

$$\frac{d{I}_i{\omega }_i}{dt}=\sum _{np}{T}_c+\sum _{nb}{T}_c$$

(2)

where $m_i$, $v_i$, $I_i$ and ${\omega }_i$ are the mass, velocity, mass moment of inertia and angular velocity of particle i, respectively. $F_g$ is the gravitational force acting on particle i. $F_c$ and $T_c$ are, respectively, the contact force and the contact torque acting on particle i due to its neighbor particles $np$ or neighbor boundaries $nb$ that are in contact with particle i. The contact forces and torques are calculated using a spring-dashpot model. For these simulations, the non-linear Hertz–Mindlin model was used, where the forces and torques are calculated as a function of the parameters stated in

Table 1. Used parameters for the discrete element method (DEM) simulation.

Property	Unit	Value
Diameter	m	0.025
Density	kg/${\mathrm{m}}^3$	2500
Coefficient of restitution, normal	−	0.9
Coefficient of restitution, tangential	−	0.5
Poisson coefficient	−	0.0235
Young’s modulus	Pa	78.5E9
Sliding friction coefficient	−	0.2
Rolling friction coefficient	−	0.002

. Further and more detailed information about the modeling can be found in the publication by Cundall and Strack [17] and in the user documentation of Simcenter STAR-CCM+ [18].

The DEM simulation was stopped when all particles were in rest. This was assumed as soon as the velocity for each particle dropped below $v_i={10}^{-6}$ m/s.

2.2. Mesh Generation

The mesh generation particle data (diameter and position of the centroid of each sphere) of the converged DEM simulation are transferred into a CAD geometry. To avoid the contact point problem, the local flattening method is used [15]. Figure 1 shows a small section of the generated mesh in the fluid domain; the mesh of the solid region is not displayed for the sake of clarity. The resulting gap width between two particles is equal to one percent of the particle diameter.

The meshing parameters are given in

Table 2. Settings of the meshing parameters.

Property	Value
Particle Diameter $d_p$	0.025 m
Base size (bs)	d
Surface edge size (ses) on spheres	0.04–0.10 $d_p$
Surface grow rate	1.3
Number of prism layers	2
Thickness of the prism layers	0.03 $d_p$
Minimum distance between two surfaces	0.12 ses

. These values are taken from our previous investigations [10,15] where the resulting mesh has proven to give accurate and mesh independent results. With these mesh settings the surface of one spherical particle is discretized with approximately 1000 surface elements. A typical mesh for a $D/d_p$-ratio = 6.8 consists of approx. 10 million fluid cells and 6 million solid cells. To minimize the influence of the inlet and outlet boundary condition, the mesh is extended by extruding the inlet by 8 and the outlet by 15 particle diameters. A schematic reactor overview is given in Figure 2.

2.3. CFD

For the CFD simulation, the finite volume method is used. The gas phase is treated as an ideal gas with constant physical properties. Constant properties are also assumed for the solid phase and are given in

Table 3. Parameters and boundary conditions used for the computational fluid dynamics (CFD) simulations.

Geometrical Properties
Particle diameter	0.025 m
Bed height	0.55 m
Diameter ratio	2.8, 4.8, 6.8, 8.8
Fluid Properties (Air)
Dynamic viscosity	1.85508E-5 Pa s
Molecular weight	28.9664 g/mol
Specific heat	1003.62 J/kg K
Thermal conductivity	0.0260305 W/m s
Solid Properties (Glass)
Density	2500.0 kg/${\mathrm{m}}^3$
Specific heat	750.0 J/kg K
Thermal conductivity	1.4 W/m s
Boundary Conditions
Pressure outlet	0.0 Pa
Inlet temperature	298.0 K
Heated wall temperature	373.0 K
Inlet velocity	0.01 m/s, 0.1 m/s, 1.0 m/s, 10.0 m/s
Particle Reynolds number	16, 160, 1600, 16,000

. The simulations were performed for a variety of tube to particle diameter ratios and particle-based Reynolds numbers $Re=\left(u_0{d}_p\rho \right)/\mu$. When applicable, turbulence was modeled using the realizable k-$\epsilon$-model with “all $y^{+}$ wall treatment”, which means that depending on the local $y^{+}$ value either a low-Reynolds formulation ($y^{+}<1$), a wall function ($y^{+}>30$) or a blend of both models ($1<y^{+}<30$) is used. For all simulations where heat transport is considered, a constant temperature of $T=100{ }^{°}$C is applied to the tube wall in the region of the packing, marked as ‘heated area’ in Figure 2. The particle conductivity is taken into account, solving the energy equation for the fluid and the solid phase. Radiation effects are neglected so that the following set of steady-state equations are solved.

Conservation of mass:

$$\nabla \left(\rho v\right)=0$$

(3)

with density $\rho$ and velocity v.

Conservation of momentum:

$$\nabla \left(\rho vv\right)=-\nabla pI+\nabla (T+T_t)$$

(4)

with unit tensor I, pressure p and the viscous and turbulent stress tensor T and $T_t$ respectively, written as:

$$T=-\frac{2}{3}\mu \nabla vI+2\mu D$$

(5)

and

$$T_t=-\mu \nabla v^{\prime }+\frac{2}{3}\rho kI$$

(6)

with dynamic viscosity $\mu$, fluctuation velocity $v^{\prime }$, turbulent kinetic energy k and the deformation tensor D:

$$D=\frac{1}{2}\left(\nabla v+{\left(\nabla v\right)}^T\right)$$

(7)

Conservation of energy:

$$\nabla \left(\rho hv\right)=-\nabla \left(pI\right)+\nabla \left(Tv\right)+\nabla \left(k\mathsf{\Theta }\right)$$

(8)

with enthalpy h, conductivity k and temperature $\mathsf{\Theta }$. Obviously for the solid phase Equations (3) and (4) are not solved, and Equation (8) reduces to

$$0=\nabla \left(k_{solid}\mathsf{\Theta }\right)$$

(9)

with the solid conductivity $k_{solid}$. All physical properties and boundary conditions used in the simulations are given in Table 3.

2.4. Presentation of the Investigated Wall Structures

To understand the effect of macroscopic wall structures on the properties of the fixed bed, six different wall structures and one tube without wall structure were investigated. These structures are depicted in Figure 3 and consist of spherical caps derived from spheres with the same particle diameter $d_p=2.5$ cm as the particles in the packing. R7 is a plain tube without any wall structure and an ideal smooth wall is assumed. R4–R6 are tubes with wall structures with a regular arrangement of the caps like those used by Zobel et al. [10] in a previous study. The wall structure and the fixed bed particle arrangement of R1 is derived from the inner zone of a fixed bed with a large tube to particle diameter ratio ($D/d_P\approx 20$) where the radial void fraction profile, averaged in circumferential direction and over the total height of the packing, already reached a constant value. This inner zone has been cut out by a cylinder. This means that all particles lying completely inside the cylindrical cut were kept and particles, which were cut by the cylinder wall were used as wall structure. The result is an absolutely irregular wall structure. For a better understanding, the procedure is outlined in Figure 4. R2 is based on the wall structure of R1, but the fixed bed is newly generated with an additional DEM simulation. R3 finally is based on the wall structure of R1 but all spheres whose centroid is located on the inner side of the tube are removed.

2.5. Post-Processing

This section shortly describes the procedures used to extract relevant data out of the CFD simulation.

The bed voidage $\epsilon$ is determined by summing up the volume $V_{f,i}$ of all N fluid cells in the geometrical range of the fixed bed divided by the volume of an empty tube with the same length L and diameter D. To minimize entry and outlet effects, the first and last three particles layers are ignored.

$$\epsilon =\frac{{\displaystyle \sum _i^N}{V}_{f,i}}{\frac{\pi }{4}{D}^{2}L}$$

(10)

The axial and radial void fraction distribution is calculated by using an area-based approach. For this, plane cross sections normal to the z-axis and cylindrical sections in the radial direction are created with an offset of 1mm respectively, leading to:

$$\epsilon \left(z\right)=\frac{A_f\left(z\right)}{\frac{\pi }{4}{D}^2}$$

(11)

and

$$\epsilon \left(r\right)=\frac{A_f\left(r\right)}{2\pi rL}$$

(12)

where $A_f\left(z\right)$ is the area of the fluid phase evaluated on cross sections along the tube axis while $A_f\left(r\right)$ is evaluated on cylindrical sections. The same cross and cylindrical sections are used to determine the radial and axial velocity and temperature profiles.

To estimate how the wall structure affects the thermal performance of the tube, the overall heat transfer coefficient U is determined with the total wall heat transfer $\dot{Q}$, heated tube surface area A, the logarithmic temperature difference $\mathsf{\Delta }T$, which is calculated with the inlet temperature $T_{inlet}$, the mass flow averaged outlet temperature $T_{out}$ and the heated wall temperature $T_{wall}$:

$$\dot{Q}=UA\mathsf{\Delta }T=UA\frac{T_{inlet}-T_{out}}{ln\left(\frac{T_{wall}-T_{inlet}}{T_{wall}-T_{out}}\right)}$$

(13)

2.6. Effect of Contact Point Treatment on Heat Transfer

Packing generation with DEM and the local flattening meshing strategy have already been validated regarding void fraction distribution, pressure drop and velocity field [15,19].

Slavin et al. [20,21] stated that the thermal conduction through the particle–particle and the particle–wall contact is negligible as long as the pressure is low, the bed itself is uncompressed and the particles are not too plastic. All three conditions are fulfilled for our setup so it can be assumed that the influence of the meshing induced local gaps on the calculated heat transfer is negligible.

To validate this hypothesis, the studies by Nijemeisland [22] (later republished in [23]) are used. The experimental setup consists of a single tube ($D=50.8$ mm) filled with 44 regularly arranged monodispersed nylon spheres ($d=25.4$ mm) and a steam-heated wall as shown in Figure 5. An air flow enters the tube and in a first unheated zone of the packed bed the flow profile develops while in the second zone the air is heated through the wall with a constant temperature ($T=100{ }^{°}$C). Temperature is measured with a thermocouple cross at the outlet of the packed bed. To cover different flow regimes, the inlet velocity was varied, resulting in particle Reynolds numbers in the range of $Re=373$ up to $Re=1922$. More information about the experimental setup can be found in the original publication.

For the simulation, three different geometries were generated that differ regarding the contact point treatment. The first geometry makes use of the global shrinking method, where each particle diameter is reduced by a certain amount resulting in a gap at the former contact points. The same shrinking percentage of 1% as in Nijemeisland and Dixon [23] is used, and therefore it is referred to as ‘Nijemeisland’. The second geometry is based on the bridging method with a bridging tube diameter of 10% of the sphere diameter d and is referred to as ‘Ookawara’. The third geometry is based on the local flattening method according to the description above and is referred to as ‘Eppinger’. To estimate the effect of the boundary layer treatment, the mesh resolution is varied between two layers and ‘all $y^{+}$ treatment’ and seven layers and a ‘Low Reynolds model’. The resulting polyhedral mesh contains approximately 700,000 cells for the model with two layers and 1.5 million cells for the model with seven layers. The physical properties are set according to

Table 4. Settings for the heat transfer validation simulation in Section 2.6.

Property	Value
${\lambda }_{fluid}$	0.026 W/(m·K)
${\lambda }_{solid}$	0.242 W/(m·K)
${\rho }_{fluid}$	1.225 kg /${\mathrm{m}}^3$
${\rho }_{solid}$	1300 kg /${\mathrm{m}}^3$
$c_{p,fluid}$	1003.6 J/(kg·K)
$c_{p,solid}$	1000 J/(kg·K)
$T_{inlet}$	298 K
$T_{wall}$	383 K

.

In Figure 6 detailed results of the dimensionless temperature $\mathsf{\Theta }=(T-T_{inlet})/(T_{wall}-T_{inlet})$ for Re = 986 are shown.

Each dot in Figure 6a–c represents a temperature value approximately 5 mm above the top layer. The numerical results show for all three methods a reasonably good agreement with a slight spread around the experimental range of data. The contact point treatment according to Nijemeisland and Ookawara slightly underpredicts the temperature profile, especially towards the center of the tube. This is not observed when the local flattening method is used. Similar results are found for Re = 373, 1724 and 1922. For an easier comparison of the different simulation results a least square best-fit curve for $\mathsf{\Theta }$ is used with

$$\mathsf{\Theta }=\frac{A}{log\left(\frac{B}{r/R}\right)}+C$$

(14)

and fitting parameters A, B and C. In Figure 6d it can be seen that all models predict the temperature profile in the near wall region quit accurately. Towards the center of the tube all models underpredict the temperature with the largest deviation for the global shrinking method while the local flattening with two layers gives the most accurate result. Similar results are also found for Re = 373, 1724 and 1922, respectively, which are depicted in Figure 7a–c:

The global shrinking method always underpredicts the temperature profile, while the local flattening method gives a slight overprediction for high Re numbers. For high Re numbers the bridging method gives the most accurate results although the particles are bridged with a cylinder with a relatively large diameter, which increases the particle heat conduction. It can be concluded that for high Re numbers the heat transport by conduction does not play a significant role, at least in this case, where the conductivity of the solid phase (nylon) is low. With respect to the boundary layer resolution the different contact point treatments are leading to different conclusions. The accuracy of the global shrinking method increases with seven layers although the temperature profile is still under-predicted. For the two other methods the two-layer approach gives better results. One reason for this might be the variation of $y^{+}$-values at the spheres surfaces that occurs due to the constant layer (cell) thickness: While it is known that the ‘all $y^{+}$ treatment’ for two-layers is able to capture $y^{+}$ values below 1 and above 30 accurately, the ‘low Re’ models require $y^{+}$ values below one, which cannot easily be guaranteed and requires a solution-adapted near-wall mesh.

As an outcome of this pre-study, we decided to use the two-layer model with local flattening because:

• the local flattening allows an automated and efficient workflow to study different configurations and
• the two-layer approach gives reasonable accurate results while keeping the mesh count and, therefore, the runtime, low.

3. Results

3.1. Preliminary Investigation

Firstly, a preliminary study was conducted to identify the most promising wall structures in terms of homogeneity of the radial void fraction distribution. The pre-study was done for all structures shown in Figure 3 with a $D/d_p$-ratio of 6.8.

In Figure 8 the circumferential and over the bed height averaged radial void fraction distribution is shown. As expected for R1 the void fraction is constant over the radius while for R7, the tube without any wall structure, the well-known oscillating void fraction profile is obtained. The profile for R2 (same wall structure as R1, but particles are refilled) is also pretty homogeneous and the typical maxima and minima of the radial void fraction profile cannot be detected. Furthermore, the bed voidage is slightly higher than for R1 since during the filling process the imposed wall structure prevents the same dense packing as for R1 and additional voids are occurring. The profiles for R3–R6 are similar to the one of R7 but the minima and maxima for the curves are not so pronounced and the void fraction at the wall is lower, around $\epsilon =0.7$. Comparing these four wall structures, it is obvious that the small caps of R6 have the lowest impact on the void fraction distribution while comparing R3–R5 with each other no significant difference can be found.

To address the question how a more homogeneous radial void fraction distribution affects the radial transport in the gas phase, a comparison of the local radial velocity in a plain section normal to the tube axis is shown in Figure 9. The range of the legend is adjusted for a better comparison in such a way that the radial velocity in red zones is equal or higher than $u_r=0.2$ m/s and in blue zones equal or lower than $u_r=-0.2$ m/s. It can be clearly seen that, with wall structures, the zones where the radial velocity exceeds the maximum and minimum is much larger. As a result, it can be expected that the lateral mixing is increased due to the higher convective transport.

To give a quantitative comparison for the different wall structures the normalized absolute radial velocity over the radial wall distance is shown in Figure 10. To remove the effect of different bed porosities ${\epsilon }_b$, normalization is done with the interstitial velocity $u_{\epsilon }=u_0/{\epsilon }_b$. All radial velocity profiles start at the wall with zero value due to the no-slip condition and show a steep increase. For R1, the maximum value is reached after approximately half a sphere diameter away from the wall and then stays constant. For reactor R7, the radial velocity magnitude oscillates around a mean value, which is significantly lower (around 30%) than that for R1. Additionally, the oscillation corresponds with the oscillation of the radial void fraction with a slight shift in the direction towards the wall. These oscillations are also found for R3–R6, but on a slightly higher level with respect to the mean averaged radial velocity magnitude (around 15% higher). The refilled configuration R2 shows a similar mean value but due to the more homogeneous void fraction distribution no distinct maxima and minima can be detected.

As a result of this preliminary study it could be shown that with wall structures the radial void fraction profile can be homogenized and as a result the radial velocity is increased which increases the lateral mixing. To quantify these effects for different flow conditions we focus in the following on only three of these wall structures. Configuration R7 without wall structures is used as a lower limiting case, R2 serves as an upper limiting case. This configuration and not R1 is chosen because R2 can be manufactured in contrast to R1, although the effort would probably be high. Finally because the wall structures R3–R6 show quite similar behavior, R5 is selected to represent these structures. This wave like wall structure was also used by Zobel et al. [10] and can be manufactured relatively easily.

In the subsequent sections these three wall structures are compared for four different tube to particle diameter ratios $D/d_p$ = 2.8, 4.8, 6.8 and 8.8 and for different particle Reynolds numbers $R{e}_p$ = 16, 160, 1600 and 16,000. For $D/d_p$ = 2.8 it was impossible to refill R2 because of jamming of the particles caused by the wall structure in combination with the narrow tube diameter; therefore, no data are presented.

Furthermore, the name for the different configurations are changed for the sake of clarity and readability. R7 without any wall structure is named plain wall reactor (PWR), R5 with the oscillating wall structure is named oscillating wallstructure reactor (OWR) and R2 with the wall structure derived from an ideal void fraction distribution is named ideal wallstructure reactor (IWR).

3.2. Bed Voidage

An overview over the bed voidage is listed in

Table 5. Bed void fraction.

	Diameter Ratio
	2.8	4.8	6.8	8.8
Plain Wall Reactor (PWR)	0.524	0.451	0.443	0.421
Oscillating Wall Reactor (OWR)	0.528	0.464	0.442	0.423
Ideal Wall Reactor (IWR)	–	0.445	0.435	0.421

. All configurations show a decreasing bed voidage for an increasing tube to particle diameter ration. This effect is reported by several authors like Dixon [24] and is caused by the reduced impact of the confining wall on the packing structure. A significant impact of the different wall structures on the bed voidage can only be observed for a diameter ratio of $D/d_P=4.8$. Compared to PWR, the OWR shows a slightly higher void fraction (+3%), while for the IWR, it is slightly lower (−3%).

The bed voidage gives no information about the spatial distribution of the particles. Therefore, the axial and radial void fraction distributions for the different wall structures and diameter ratios are displayed in Figure 11 and Figure 12. For a diameter ratio of $D/d_P=8.8$ the axial void fraction distribution shows the well-known oscillating trend starting with a value of $\epsilon =1$ and ending in a more or less constant value after a distance of about three to four particle diameter away from the bed entry. With decreasing diameter ratio the oscillation around the average value increases indicating that the smaller tube diameter determines more and more the particle position and prevents the generation of a fully random packing. Especially for $D/d_P=2.8$ a quite regular stacked bed of particles is found, which ends in an overall oscillating void fraction distribution. A significant impact of the wall structures on the axial void fraction distribution cannot be observed for all investigated cases.

The radial void fraction distribution depicted in Figure 12 shows a significant impact of the different wall structures on the particle arrangement for all diameter ratios except for $D/d_P=2.8$. Here, similar to the axial void fraction distribution the confining wall of the small tube determines the position of the particles. For larger $D/d_P$-ratios the radial void fraction profile shows the same trend as discussed in the preliminary study section: Compared to PWR the OWR shows a more homogeneous profile with less pronounced peaks, while the positions of the minima and maxima are nearly unaffected. The irregular wall structure of the IWR leads to an almost constant radial void fraction distribution. The void fraction value at the wall decreases from $\epsilon =1$ for PWR over $\epsilon \approx 0.6$ for OWR to $\epsilon \approx 0.4$ for IWR for all $D/d_P$-ratios.

3.3. Velocity

To understand the impact of the homogenization of the radial void fraction distribution on the velocity field Figure 13 shows the axial progression of the cross-section averaged magnitude of the radial velocity

$$\left|u_r\right|=\frac{1}{{\displaystyle \sum _i^N}{A}_{f,i}}\sum _i^N\left|u_{r,i}\right|·A_{f,i},$$

(15)

normalized with the interstitial velocity magnitude

$$u_{\epsilon }=\frac{1}{{\displaystyle \sum _i^N}{V}_{f,i}}\sum _i^N{u}_i·V_{f,i}.$$

(16)

$|u_r|$ and $u_{\epsilon }$ were evaluated on cross sections normal to the z-axis along the whole reactor. The fixed bed is situated at $0\le z/L\le 1$ and $|u_r|/u_{\epsilon }$ can be interpreted as the intensity of the convective radial mixing.

Figure 13 shows that for all investigated Reynolds numbers an increasing convective radial mixing can be observed for OWR and IWR compared to PWR. Only for a diameter ratio of $D/d_P=2.8$ this effect does not exist because of the same reason as discussed earlier. To get more quantitative information

Table 6. Axially averaged mean radial velocity magnitude.

	Diameter Ratio
	2.8			4.8			6.8			8.8
Re	PWR	OWR	IWR	PWR	OWR	IWR	PWR	OWR	IWR	PWR	OWR	IWR
16	0.181	0.177	–	0.210	0.246	0.297	0.226	0.274	0.295	0.240	0.290	0.298
160	0.202	0.195	–	0.237	0.261	0.306	0.242	0.283	0.302	0.256	0.296	0.303
1600	0.205	0.197	–	0.230	0.248	0.281	0.219	0.272	0.278	0.233	0.257	0.262
16,000	0.188	0.183	–	0.204	0.223	0.257	0.196	0.209	0.256	0.210	0.242	0.248

shows the axially average values of $|u_r|/u_{\epsilon }$ in the range of $0\le z/L\le 1$.

Comparing the effect of different Reynolds numbers for each configuration individually a maximum of the convective radial mixing around $Re=160$ is found. For higher and lower values the lateral mixing is decreasing. Additionally, for most Re the radial mixing increases with rising diameter ratio for PWR and OWR. This can be explained by a reduced significance of the bypass effect and, therefore, also of the wall structure closed to the confining wall when the diameter ratio is raised. For the IWR this effect is not distinct. This is an evidence that the ideal wall structure leads to a more or less homogeneous radial void fraction distribution which prevents bypass effects. Table 6 illustrates also that there is a significant increase in the lateral mixing when a wall structure is used. Compared to PWR the OWR shows an increase of radial mixing of 6–21%. The benefit of the IWR lies in the range of 12–41%.

3.4. Heat Transfer

Figure 14 shows the temperature distribution in a cross section normal to the x-axis. The four plots on the left Figure 14a–d show the well-known parabolic temperature profile for tube reactors with varying Reynolds Number, although the ideal parabolic profile is disturbed by the presence of the solid particles. It can clearly be seen that with increasing velocity the core temperature along the z-axis stays longer with the inlet temperature. For the lowest velocity (Figure 14a) immediately at the packing entry the temperature rises. In contrast, for the highest investigated velocity (Figure 14d) even at the end of the bed the core temperature has not drastically changed. Comparing the three different configurations at $Re=1600$ in Figure 14e–g it can be clearly seen that the overall temperature for the plain wall configurations is the lowest. OWR shows a higher overall temperature, especially in the near-wall region. This is even more pronounced for IWR and leads to volumes close to the wall where the wall temperature is almost reached. This indicates a kind of dead zone or zone with a significantly reduced convective heat transport.

To determine the impact of the improved convective radial mixing on the heat transfer the thermal entry length has been estimated. Therefore, the axial core temperature profile in the range of the fixed bed has been extracted using line probes. After that the thermal entry length was calculated searching for the first axial coordinate where the dimensionless temperature $\mathsf{\Theta }$ reached a value of $\mathsf{\Theta }\ge 0.01$. The calculated thermal entry lengths are listed in

Table 7. Calculated thermal entry length.

	Diameter Ratio
	2.8			4.8			6.8			8.8
Re	PWR	OWR	IWR	PWR	OWR	IWR	PWR	OWR	IWR	PWR	OWR	IWR
16	0.006	0.000	–	0.000	0.006	0.000	0.000	0.000	0.000	0.000	0.000	0.000
160	0.044	0.006	–	0.033	0.033	0.011	0.105	0.078	0.066	0.170	0.159	0.165
1600	0.050	0.022	–	0.116	0.116	0.072	0.281	0.221	0.122	0.461	0.428	0.352
16,000	0.189	0.105	–	0.237	0.199	0.116	0.447	0.354	0.304	$>L$	$>L$	$>L$

.

The results show that firstly no thermal entry length can be determined for the lowest Reynolds number (Re = 16). In this case the conductivity driven heat transfer becomes dominant. Secondly the thermal entry length exceeds the bed length for Re = 16,000 and $D/d_P=8.8$. The remaining data in Table 7 show a trend of increasing thermal entry lengths with rising diameter ratio and increasing Reynolds number as expected. Concerning the impact of the wall structures on the heat transfer, the data prove a significant reduction of the thermal entry length up to 66% for the IWR compared to PWR ($D/d_P=4.8$, $Re=160$). The data also show that in general the reduction of the thermal entry length is higher for IWR than for OWR indicating the importance of the homogeneous void fraction distribution. Finally, as discussed in Section 3.3, the effect of the wall structure is less pronounced for higher $D/d_P$ values because of the lower importance of the bypassing close to the wall.

To get a better understanding of the temperature field, Figure 15 displays the circumferential averaged dimensionless radial temperature profile at four different axial positions within the fixed bed. Because of the huge amount of data, this is done for different tube to particle diameter ratios at a constant Reynolds number of Re = 1600. Similar trends were found for the other Reynolds numbers, too. For all configurations and all positions in the bed the PWR shows the lowest radial temperature profile and IWR the highest while OWR is somewhere in between. For both configurations with wall structure the temperature profile is smooth while for PWR steps in the profile can be seen, especially in Figure 15a,b. Again, this can be explained by the particle arrangement caused by the confining wall and its resulting radial void fraction distribution.

3.5. Heat Transfer Coefficient

To estimate how the wall structure affects the thermal performance of the tube, the overall heat transfer coefficient U is determined with Equation (13). The product of the heat transfer coefficient U and the heated surface area A is a measure of how much heat can be transferred for a given temperature difference. Figure 16 shows UA for all investigated cases, normalized with the values of the plane tube (PWR) for a better comparison.

In Figure 16a the values are shown for $D/d_p$ = 2.8. As mentioned earlier, IWR could not be investigated due to blockage during the refilling of this narrow tube. The oscillating wall structure OWR shows an increase of UA up to 25% with increasing Reynolds number. A maximum is reached for Re = 1600 and above.

Figure 16b depicts the results for $D/d_p$ = 4.8. PWR shows similar values as for $D/d_p$ = 2.8. The tube with the ideal wall structure IWR shows over the whole range of Reynolds numbers a higher value for UA, but like for PWR a maximum is reached at Re = 1600 and a similar value is found for Re = 16,000. The value UA is more than doubled for these Reynolds numbers compared to PWR. This indicates that the effect of the wall structure is more pronounced for the turbulent flow regime and confirms that an effective disturbance of the by-passing improves the thermal performance significantly.

Similar results are found for $D/d_p$ = 6.8 in Figure 16c and for $D/d_p$ = 8.8 in Figure 16d. It should be noted that for the largest $D/d_p$-ratio the improvement for PWR as well as for IWR is slightly reduced. This can be explained by the reduced importance of the near wall region with increasing tube diameter. For larger $D/d_p$-ratios the heat transport in the center of the bed is getting more and more important and starts to dominate the performance. Since the local void fraction between the different tubes aligns to each other after approximately four particle diameters as shown in Figure 12, it can be expected that the UA values are equalized for much larger $D/d_p$-ratios.

4. Conclusions and Outlook

Regarding the results and discussion of this study, the following conclusions can be drawn:

• The local flattening meshing strategy [15] is appropriate to simulate the heat transfer in fixed bed reactors.
• It is possible to homogenize the radial void fraction distribution with the use of an appropriate macroscopic wall structure. This leads to a significantly improved radial convective mixing.
• The improvement of the lateral mixing is dependent on the Reynolds number and the tube to particle diameter ratio. There seems to be an optimal Reynolds number where $|u_r|/u_{\epsilon }$ achieves a maximum. In this study, an increase between 6 and 41% could be reached depending on the flow condition and the type of macroscopic wall structure.
• The better radial mixing leads to a shorter thermal entry length and a higher radial temperature profile within the reactor. Depending on the flow condition and on the configuration the thermal entry length is reduced between a few and 66%.
• The thermal performance can be improved by up to 120% depending on the wall structure and flow regime. The improvement is the highest for industrial relevant Reynolds numbers.

In future works the technical implementation of macroscopic wall structures with the help of 3D printing will be investigated and the promising results of this study should be validated with a reactive system where the advantages of the macroscopic wall structure on the heat transfer and as a consequence on conversion and yield can be quantified.