2.2.2. Pseudo-Homogeneous *λ*eff,r(*r*)-*u*z(*r*) 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\_{\rm f} c\_{\rm p,f} u\_{\rm x}(r) \frac{\partial T}{\partial z} = \frac{1}{r} \frac{\partial}{\partial r} \left( \lambda\_{\rm eff,x}(r) r \frac{\partial T}{\partial r} \right) + \lambda\_{\rm eff,x}(r) \frac{\partial^2 T}{\partial z^2}. \tag{21}
$$

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

$$T = T\_{\rm W} \tag{22}$$

As reviewed by Dixon [31], multiple models exist to determine *λ*eff,r(*r*). 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 *λ*eff,r(*r*) = *f*(*r*). Ahmed and Fahien [55] defined the wall region to be 2 *d*p,v thick and used a cubic dependency for *λ*eff,r(*r*) 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 *λ*eff,r in the center of the bed, at the tube wall, and at the interface of both regions. Contrary, Gunn et al. [57–59] used a constant value for *λ*eff,r in the bulk region and assumed a quadratic dependency of *T*(*r*) in the wall region. They defined the wall region to be 0.3 *d*p,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 *λ*eff,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 *λ*eff,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 *λ*eff,r(*r*) = *f ρ*f, *c*p,f, *d*p,v,*ε*0,*ε*w,*ε*(*r*), *u*0(*r*) 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\_{\rm eff,r}(r) = \lambda\_{\rm eff,r}^0 + \frac{u\_{0,\rm c} d\_{\rm P,v} \rho\_{\rm f} c\_{\rm P,f}}{K\_{\rm W}} \cdot f(R - r), \tag{2.3}
$$

using:

$$f(R - r) = \begin{cases} \left(\frac{R - r}{k\_{\rm f,w}d\_{\rm p,v}}\right)^{n\_{\rm f,s}} & \text{if } 0 < R - r < k\_{\rm f,w}d\_{\rm p,v} \\ 1 & \text{if } k\_{\rm f,w}d\_{\rm p,v} < R - r < R \end{cases} \tag{24}$$

was used as the basis to determine *λ*eff,r(*r*). In Equation (23), *K*<sup>w</sup> 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*f,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*f,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* = [0.1 : 0.1 : 1.1] were extracted from the particle-resolved CFD simulations for the simulation with *Re*<sup>p</sup> ≥ 500. For the lowest investigated Reynolds number of *Re*<sup>p</sup> = 100, the radial temperature gradients flattened out quickly. Therefore, in this case, only the temperature profiles in the range of *z* = [0.1 : 0.1 : 0.3] 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 (*Re*<sup>p</sup> ≥ 500) and 300 (*Re*<sup>p</sup> = 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.

**Figure 3.** Visualization of the geometry and the mesh for a 6-hole tablet (**left**) and a 10-hole cylinder (**right**).

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 Nm3/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: *<sup>λ</sup>*s[*W*/(*m K*)] = 8.195 · exp 1.188 · <sup>10</sup>−<sup>3</sup> · *<sup>T</sup>* . 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 Θ = (*T* − *T*0)/(*T*ref − *T*0), whereas *T*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* ≈ 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* ≤ 0.35. For a flow rate of 15 N m3/h, also deep in the bed, a good agreement with the experimental data was found. However, at higher flow rates, for *z*/*h* ≈ 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.

**Figure 4.** Axial profiles of dimensionless temperature at different flow rates for a 6-hole tablet (**left**) and a 10-hole cylinder (**right**). Comparison of CFD results (scattered cloud) against experimental data (big symbols).

#### **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* = *Q*˙ w/ *A*wΔlog*T* , using Δlog*T* = (*T*out − *T*in)/ log((*T*<sup>w</sup> − *T*in)/(*T*<sup>w</sup> − *T*out)), was used to compare the different designs, whereby *U* was evaluated for an axial threshold of the reactor that fulfilled the criterion 0.0 ≤ Θcore ≤ 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. A schematic drawing of the setup is presented in Figure 5.

**Figure 5.** Schematic drawing of the numerical setup, including a snippet of the mesh. Exemplarily shown for the loose packing of spheres.

**Table 2.** Summary of important results and simulation parameters.

