Particle-resolved computational modeling of hydrogen-based direct reduction of iron ore pellets in a fixed bed. Part II: Influence of the pellet sizes and shapes

Full-fledged computational modeling of Direct Reduction reactors encompasses single-pellets models and the step-wise scaling up to industrial-scale reactors. This study delves into the particle-resolved computational modeling of hydrogen-based direct reduction of 0.5 kg iron ore fixed-beds, as a scale-bridging step and focusing on the influence of pellet sizes and shapes. To investigate the effect of particle sizes, two beds with particles sized 10.0 – 12.5 and 12.5 – 16.0 mm were characterized and numerically reconstructed using the discrete element method. While to assess the effect of particle shapes, a third numerical bed was generated directly upon a high-resolution CT-scan of a real bed. Through reactive CFD simulations, we investigated the reduction process in these three beds using hydrogen as the reducing gas. Our findings reveal that the bed structure significantly impacts the reduction efficiency and overall conversion degree. This study emphasizes the importance of accurate bed reconstruction and provides critical insights for optimizing the hydrogen-based direct reduction process in industrial applications.


Introduction
The steelmaking industry is mostly dominated by the blast furnace/ basic oxygen furnace (BF/BOF).However, the inherent environmental and economic drawbacks of the process prompted the development of alternative processes such as mini blast furnaces, smelting reduction, and direct reduction (DR) processes [1].Among the alternative processes, the DR process with green hydrogen is considered to be the most promising technology for decarbonizing the steel industry [2][3][4].In 2022, the global direct reduced iron (DRI) production was 127.36 million tons (Mt), and considering the last six years, the DRI output has grown by almost 55 Mt, that means around 75% [5].
Over the past century, many different direct reduction processes have been developed, but they have been abandoned due to economic and technical reasons [6].To date, there are two popular groups of DR processes in operation, namely coal-based and gas-based processes.In 2022, 28% of DRI productions are from coal-based rotary kiln processes [5], whereas over 70% of production comes from gas-based shaft furnaces, and a small fraction of gas-based DRI is produced using fluidized bed reactors.The shaft furnace process is dominated by the ENERGIRON and Midrex processes, accounting for around 12.1% and 57.8% of total DRI production in 2022, respectively [5].
The iron ore concentrates, obtained through the enrichment processes involving size reduction, need to undergo agglomeration in a pelletizing plant [7].As efficient reduction of iron ore is highly favored by higher reducing gas permeability and a high amount of solid-gas interactions, the iron ore raw material requires certain physical characteristics.The general requirements are to be of specific size, metallurgically reducible, have no undesirable elements, have a non-sticking tendency, have an acceptable chemical composition, etc.To ensure better gas permeability, the size of the pellets needs to have a proper distribution [6].The shape of industrial iron ore pellets used in the DR process varies based on several factors, and it plays a crucial role in the efficiency and performance of the reduction process [7].Spherical or near-spherical pellets are ideally expected in the DR process to promote uniform heat distribution, better gas flowability, and to reduce the risk of clogging or uneven distribution in various stages of steel production.However, iron ore pellets for DR are industrially produced from natural hematite grains (irregular, ~20 μm) with a generally rough spherical form [8]. Typical pellets are 9-16 mm in diameter with a varying degree of sphericity [9].Achieving near-perfect spheres is impractical from the actual pelletization process.Pelletization involves agglomerating fine iron ore particles into larger, cohesive pellets, typically using binders.The type of equipment and technology used for pelletization can affect the final shape of the pellet [7].Other factors include binder and additives such as biomass [10], moisture content during pelletization, and the drying process, among other factors.
In fixed-bed reactors where the tube-to-particle diameter ratios λ = ( D /d p ) is low (≤15), the irregular bed structure plays a crucial role in influencing fluid dynamics, thereby affecting heat and mass transfer [11].Rodrigues et al. [12] investigated how the particle shape can affect the pore morphology of a randomly generated packed bed.The group proposed a relationship between bed tortuosity and particle sphericity.The bed tortuosity increases parabolically as the sphericity of the particles decreases from 1 to 0.8; after that, the tortuosity is constant.It is worth noting that the tortuosity of the bed characterizes the mass transport in a porous bed and, hence, can influence the overall gas-solid reaction speed.The authors also investigated the impact of the sphericity factor on the radial distribution of the bed porosity, and they concluded that the bed with higher spherical particles exhibits strong oscillations of porosity along the radius, whereas lower sphericities result in irregular and faded oscillations.
Phase-averaged CFD models may be inadequate for representing particle-level interactions and intricate fluid dynamics, especially when considering local flow effects [11].This is especially true for tube-to-particle diameter ratios between 4 and 7 [13].This phase-averaged CFD model considers the average value of the bed porosity and makes no clear distinction between the phases.The coupled CFD-DEM model combines discrete element methods for the solid phase with CFD for the fluid phase, suitable for various iron ore reduction cases [14][15][16] but does not fully resolve flow around particles.This CFD-DEM approach is also called unresolved approached [17], as the grid size used for CFD is much larger than the particles.Hence, to investigate the physical phenomena in a fixed bed reactor, the application of computational fluid dynamics (CFD) in a particle-resolved bed structure appears to be more appropriate.
The particle-resolved bed structure considers the actual arrangement of the pellets in the bed, which means kinetics and transport phenomena can be considered in the interstitial region of the pellets [11].For CFD simulations of a particle-resolved fixed bed, the initial step is to develop a geometric representation of the bed using either scanning techniques (reconstructive methods) or numerical methods.Reconstructive methods use advanced scanning technologies such as MRI and X-ray microtomography (XMT) to capture the exact shapes and positions of particles.Many works [18][19][20], including the present one, utilized this technique for geometric bed reconstruction.While this approach offers high fidelity, it is often time-consuming and computationally intensive.The numerical techniques fall into three main categories: idealized particle arrangements, random particle arrangements, and Monte Carlo methods [11].Idealized particle arrangements involve placing particles in simple geometric patterns like simple cubic (SC), face-centered cubic (FCC), and body-centered cubic (BCC).Authors like Ferng et al. [21], Lee et al. [22] applied this method to investigate turbulence and corresponding heat transfer in a bed.These methods are computationally efficient but less representative of actual conditions, especially for non-spherical particles.Random particle arrangements employ techniques like the Discrete Element Method (DEM) to simulate realistic packing by modeling particle interactions and dynamics, achieving a balance between realism and computational load.This work utilizes this random particle arrangement as second technique, beside the reconstructive method, to generate the iron ore pellet bed.Monte Carlo methods generate packings by iteratively adjusting particle positions to reduce voids.While flexible, they may sometimes produce non-physical configurations as the method does not consider particle collisions, particularly with complex particle shapes or reactor internals.To address this non-physical behavior, Caulkin et al. [23,24] extended this method by introducing a hybrid approach known as Collision Guided Packing (DigiCGP).Although this approach accounts for particle collisions, it still encounters difficulties in accurately capturing the local particle orientation, which may impact the precision of bed morphology simulations.In the first part of this two-part article [25], we have demonstrated the use of the discrete element method [26] to generate the bed and performed CFD simulation of the fixed bed reduction.In this second part here, we explore further the influence of the pellet sizes and shapes on the reduction process.This study will address the existing knowledge gaps regarding the impact of bed structure, pellet size, and shape distribution on iron ore reduction kinetics and gas flow dynamics, which require further exploration.The objectives of this work are as follows.
1. Characterization of two iron ore beds with similar weight for pellet size distribution using analytical (image analysis) and experimental method (water displacement method), determining the bed structures (axial averaged bed porosity).2. Reconstruction of those beds using a.A numerical technique (DEM method) and b.An advanced tomographic method (CT-scan) 3. Generation of 3D mesh structures based on the positions and sizes of the pellets.Another mesh will be created using a 3D CT-scan image.4. Performing CFD simulation of those reconstructed beds for the direct reduction of iron ore process using H 2 gas. 5. Conducting a comparative analysis to demonstrate how variations in bed morphology (size/shape distribution) influence the reduction process in fixed-bed reactor.Potential improvements and optimization to increase the overall reduction efficiency are also discussed.
To this attempt, two batches (noted T1 and T2), containing industrial pellets with two distinct distributions and weighing approximately 0.582 kg each, have been sampled from larger bulk supplies.The particle size distribution, shapes, and resulting fixed bed structures have been characterized using various methods, including CT-scan from the T1 batch.In addition, numerical reconstructions of the beds have been performed (T1 DEM , T2 DEM and T1 CT ) followed by CFD simulations of the reactive systems.Numerical beds denoted T1 DEM and T2 DEM follow the methodology described in the first part, where DEM simulations are conducted to settle the pellets knowing the particle size distribution as input.Then, a mesh generator is used to create an appropriate 3D mesh based on the positions and sizes of the pellets in the modeled bed, resolving skewed cells in narrow contact regions among pellets through appropriate treatments and refinements.Because this technique is limited to using spherical particles, a mesh suitable for CFD simulation -T1 CT -was also generated directly from the CT-scan, thus allowing to work on the real pellet shapes.For all three beds, the same CFD simulation framework has been utilized.It relies on an in-house transient solver developed in the open-source CFD software OpenFOAM [27] for the gaseous reduction of iron ore pellets.The solver considers reactive gas-solid interactions, gas species transport within the porous pellets, and porosity variations.Kinetic and transport data from a previously established 1D model for the direct reduction of iron ore pellets in H 2 atmospheres are applied [28,29].A comparative analysis of three different bed structures have been performed to demonstrate how different bed structures influence the reduction process of iron ore pellets in a fixed bed setup.

Materials
The iron ore pellets used in this work are commercial pellets collected from a source produced in Türkiye.Fig. 1 presents a couple of iron ore pellet samples used in this work, which show that the pellets are quite irregular in shape.
The basic characteristics, along with the major and trace elements present in the pellets, are listed in Table 1.Two different pellet bulk supplies, S1 and S2, of approximately 20 kg each were provided.The size distribution mentioned in the table has been given by the manufacturer.

Experimental construction of fixed beds
To ensure a representative sampling, a total of 0.582 kg of pellets, comprising 175 (T1) and 84 (T2) individual pellets, were randomly collected from the bulk supplies.A 72-mm-diameter Plexiglas cylinder was used to accommodate the pellets, creating two approx.7-cm-high fixed beds.The Plexiglas cylinder further served to produce a CT-scan tomographic image.

Size distribution
To gauge the pellet size distribution, the two Turkish pellet samples were analyzed.Due to their irregular spheroid shape, three different methods to calculate the pellet diameters were employed.
Methods I and II: The samples were subjected to image analysis.The pellets were spread out on a white background and photographed with an overhead camera and light.The images have been further converted to binary images in ImageJ [30].The particle analyzer tool was then utilized to extract the area, circularity, and ellipse fit of each pellet, where circularity is defined as: 4π⋅area perimeter 2 .The size distributions of both pellet samples were constructed based on the diameter from a circular fit of the pellet area and the minor axis of the ellipse fits, respectively.
Method III: A water displacement method were employed [31].The pellets of the two samples were weighted individually and then placed in a cylinder with a diameter of 7.2 cm and immersed in water until just covered.The water was strained from the pellets and weighted to determine the void volume between them, which was used to calculate the volume of the pellets themselves occupying the cylinder.From the pellet weight before the measurement and their total volume, the pellet's bulk density could be calculated.Using this density and the mass of each pellet, their individual volumes were obtained.The size distribution was then constructed by using a spherical assumption for the pellet diameters.

Average radial bed porosity
To understand the morphology of a fixed bed structure, it is common to investigate how the porosity of the bed varies across the radial direction.A non-destructive way to do that is to capture the threedimensional structure of the bed.An X-ray tomography image of the bed was created by Messtronik GmbH (St.Georgen, Germany) as sketched in Fig. 2 (a).The model was then sliced vertically through its center, traversing eight 2D planes using the Slic3r software [32], as shown in Fig. 2 (b and c).These eight planes were taken at an angle of 0 • , 20 • , 45 • , 60 • , 90 • , 110 • , 135 • , and 150 • and subsequently averaged to determine the average radial distribution of the bed porosity.

Numerical reconstruction of fixed beds
The Discrete Element Method is an explicit numerical method used for simulating the movements and interactions of particles, including collisions with walls.OpenFOAM [16] uses a soft-sphere-based approach for DEM simulation.The "icoUncoupledKinema-ticParcelFoam"-solver is used to model the interaction between particles and walls.This approach treats particles as deformable bodies capable of undergoing multiple, simultaneous collisions.Newton's law of motion is used to describe each particle's movement, considering particle-particle and particle-wall interactions as well as gravitational forces.The drag force from the gas phase (continuum phase) is neglected, resulting in the momentum conservation equation as follows: In the left part of the equation, m i is the mass of the particle i and dvi dt represents changes in the i th particle velocity over time t.On the righthand side of equation ( 1), F g and F c represent the gravitational force acting on the particles and the contact forces, respectively.F c is the sum of the particle-particle and particle-wall interactions.These contact forces are calculated using the spring-dashpot model.When the particles Fig. 1.Sample iron ore pellets.Highlighting size and shape variability.

Table 1
Physical and chemical characteristics of iron ore pellets.Two fixed beds structures with the experimental size distributions were constructed using the DEM simulation.Table 2 presents the characteristics of the two bed structures.Irrespective of the size distribution and pellet number, the total sample weight is 582 g.
The numerical parameters essential for calculating the contact forces utilized in the spring-dashpot model are presented in Table 3.The equations used to calculate particle-particle collisions are also valid for calculating particle-wall collision forces by assuming that the wall is a rigid surface.The parameters for the particle-particle and particle-wall collision models used in the DEM simulations are detailed in Table 4.Some of the coefficients used here are taken from literature and a sensitivity analysis of them helps to determine the influence of the parameters on the simulation outcomes.
Based on the particle positions and size distribution of the T1 and T2 bed structures from the DEM simulations, snappyHexMesh, an Open-FOAM utility, was used to create the mesh.An optimized volume mesh has been generated for the calculation domain, consisting of both gas and solid phases.Fig. 4 shows a detailed view of the applied mesh in an exemplary slice in the direction of the flow.A global shrinking method of 1% [33,34] was implemented to deal with the contact point problem.This method introduces a minimal distance between two contact surfaces and facilitate the mesh generator to have adequate space to create good-quality cells.This prevents cell skewness and enhances convergence.
The pressure drop across the bed is expected to be minimal, based on the setup with an inlet velocity of 0.805 m s − 1 and bed heights of 78 mm (T1 DEM , T2 DEM ), and 65 mm (T1 CT ).According to Pakov et al. [35], the pressure drop is significantly influenced by bed height, with shorter beds   like ours exhibiting minimal pressure losses.Furthermore, the study of Pashchenko et al. [36] indicates that at low inlet velocities such as 0.805 m s − 1 , the kinetic energy of the gas is insufficient to create substantial pressure drops, especially given the short bed height.Additionally, Pashchenko [37] highlights that while particle roughness can increase pressure drop, this effect is more pronounced in taller beds or at higher gas velocities.The most challenging mesh generation is for the T1 CT structure, where the CT-scanned model in section 2.2 is utilized for creating the complex internal faces.The problematic contours and holes in the CT-scanned model were examined using a mesh manipulation program called Meshmixer.These problematic regions were then repaired, making the model watertight and solid with no holes or internal baffles.The repaired 3D model was subsequently converted into an STL (stereolithography) file, which represents the surface geometry of the pellet bed.A structured background mesh has been created with the blockMesh utility in OpenFOAM.And using the snappyHexMesh, the repaired STL file has been defined as bed geometry, which allows the surface to be refined and snapped.As depicted in Fig. 4 (c), the bed region undergoes refinement at a level 1. Fig. 4 (d) shows part of the final mesh, while Fig. 5 shows the steps for producing this real bed mesh structure.To reveal the intricate surface structure, the solid and gas phase mesh regions are presented separately in Fig. 4 (d).An optimized mesh has been created for T1 CT structure with a total number of cells comparable to T1 DEM and T2 DEM beds, as listed in Table 5.
All the beds are approximately 78 mm high.There is a 10 mm empty calculation domain at the bottom of the pellet bed and a 50 mm empty calculation domain at the top of the pellet bed to encounter the backflow effects during simulation.

Reactive CFD modeling of iron-ore pellets
A native OpenFOAM transient solver [27] has been modified and applied for modeling the fixed-bed iron ore pellet reduction process.The governing equations for modeling a laminar reactive flow through the fixed-bed structure include the Navier-Stokes equations (conservation of mass and momentum) and the conservation equations of species.The equations are presented in Cartesian coordinates and are valid for the whole reactor domain including the pellets and the gas phase.The governing equations for gas-surface interactions are as follows: Mass conservation (Continuity equation) : Momentum conservation (Navier − Stokes momentum equation) : Species conservation (Species transport equation) : The section Nomenclature defines the symbols used in the governing equations mentioned above and throughout the article.More detailed discussion regarding the governing equations are done in one of our previous publications [29].The void volume fraction (α g ) is equal to 1 for the gas phase and no homogenous chemistry among the gaseous species is considered.As a result, the net mass transfer between gas and solid phases S g , the momentum sink term (M sg ) and the chemical sink terms ( ṡi ) are 0. The governing equations for the gas phase are then appears to be ∂ρ g ∂t ∂ρ g v g ∂t The effective diffusion flux j eff g,i = − ρ g D eff ∇Y g,i , where D eff is the effective diffusion coefficient.In our work, only the molecular diffusion is considered when calculating the effective diffusion, and no Knudsen diffusion is considered.This is due to the fact that the porous structure of industrial iron ore pellets is not within the nanometer range, and the DR process is carried out at atmospheric pressure [38].
The governing equations for the solid phase are given by: ∂ε ∂t The reaction-diffusion inside each pellet is solved by using a pseudocontinuum model, where each of the pellets is considered as a continuum phase and is described by the void volume fraction (α g ) and the tortuosity (τ).The transport within these porous pellets is described using the effective diffusivity.The porosity of the pellets is modeled by using the Darcy law applied to the defined cell zone, where the momentum equation considers a new source term (M sg ).This added source term in the momentum equation can then apply to flow-directed pressure drops.In this study, a reducing gas mixture containing 75% H 2 and 25% N 2 at a flow rate of 50 Nl/min at a temperature of 1173 K has been introduced from the bottom of the reactor, as depicted in Fig. 6.The boundary conditions applied for the CFD simulation are listed in Table 6.
The direct reduction of iron ore using hydrogen includes several stages with following sequence of oxides: Fe 2 O 3 (hematite) → Fe 3 O 4 (magnetite) → FeO (wüstite) → Fe (metallic iron).In order to develop a reliable chemical mechanism for the DR process, we proposed a 1D modeling framework that solves the reduction process at the pellet level, considering spherical symmetry and coupled it with a 3D CFD model [29].In this way, the concentration gradient between the outer surface of the pellet and the bulk gas is considered.Our 1D model framework has already shown maturity in predicting much more complex reduction behavior, including phenomena like cementite formation and carbon deposition [28].
The modeling framework is developed under the assumption of isothermal conditions, but in reality, the bed temperature varies with time.To consider this non-isothermal effect, we adopted an experimental temperature profile for the H 2 reduction process [39] with comparable experimental conditions.In our CFD solver, we dynamically adopted the experimental temperature at runtime to simulate the non-isothermal reduction process.The reactions parameters listed in Table 7 and the kinetic are derived in Part I of this work by validating a fixed bed experimental setup [25].For simplification purposes and ease of parameter estimation, all reactions are considered homogenous and irreversible.
The chemical source term S, used in the governing equations is calculated as follows: This is the chemical source for the first reaction in Table 7 (Fe 2 O 3 → Fe 3 O 4 ), where, C H2 is the Hydrogen gas concentration and X Fe2O3 is the mole fraction of Fe 2 O 3 .The Arrhenius relation is used to calculate the rate constants, Here, R and T are the universal gas constant and temperature, respec-  The overall conversion degree, F, is the parameter used to quantify the extent of conversion (from iron ore to reduced iron) of the bed.The overall conversion degree is defined as where the following notation is used: m 0 -initial pellet mass, m(t) mass at time t, m ∞ -theoretical mass after complete reduction.The values 0 and 1 for F indicate no conversion (Fe 2 O 3 and Gangue) and complete conversion (Fe and Gangue), respectively.

Evaluation of experimental size and shape distributions
The three methods described in section 2.2.1 are compared to establish a size distribution for the two samples (T1, T2) serving as a foundation for the DEM simulations.Fig. 7 shows the successive steps of the image analysis for method I and II.The area of each pellet is measured from a binary image (Fig. 7 (b) and (c)).Referring to a circular geometry, the areas obtained were converted into their respective diameters.
The results for both samples, T1 and T2, are visualized in Fig. 8.The distributions exhibit characteristics of normal distributions; nevertheless, they do not closely conform to the manufacturer's specified size ranges of 10.0-12.5 mm (T1) and 12.5-16.0mm (T2), with expected average sizes of 12.7 mm (T1) and 16.2 mm (T2).The images of Fig. 7 (a) and (b) show a significant non-sphericity of the pellets.Circularity, as defined above, was computed for each pellet, resulting in an average circularity of 0.81 for sample T1 and 0.84 for sample T2.Using this circularity as a correction factor for each pellet's diameter yields the distributions in green in Fig. 8, with an average size of 10.4 mm (T1) and 13.6 mm (T2).This observation suggests that the circularity could serve as a viable correction factor for pellets with substantial non-sphericity.Particularly, in cases where the manufacturer's size range is achieved by sieving pellets through various mesh sizes, the emphasis on the smallest feasible diameter becomes a primary selection criterion.The circularity factor adjustment effectively mimics this selection process for non-spherical particles.
In the second image analysis approach, this smallest feasible diameter for each pellet is selected by using the minor axis of an ellipse fit (Fig. 7 (d)).The resulting distributions are shown in Fig. 9 (blue) and agree very well within the range given by the manufacturer, with an expected mean size of 11.6 mm (T1) and 15.3 mm (T2).This is arguably a better measure, that accounts for the non-uniformity of the pellets.A third approach used the individual pellet weight and bulk density obtained from the water displacement method to calculate the volume of each pellet.Using this method, bulk densities of 3.23 g/cm 3 (T1) and 3.10 g/cm 3 (T2) were obtained.The pellet diameter is calculated using a spherical assumption, and is shown in Fig. 9 (red).The expected sizes are very similar to the elliptical fit, with an expected average size of 11.8 mm (T1) and 15.4 mm (T2).The volumetric approach confirms the overestimation of the pellet size of non-uniform pellets when only the pellet area is considered in image analysis procedures.For comparison, the extracted volumes of each approach are shown in Table 8.Using the area-diameter to generate the pellet volume with a spherical assumption largely overestimates the volume compared to the volume observed by the water displacement method.The pellet volume from the circularity corrected diameter however is greatly underestimated.
The total volume calculated using the minor axis approach closely adheres to the total volume calculated by individual volume measurement.This result is consistent with the similar size distributions (Fig. 9) calculated by those methods.Despite both results concluding that the minor axis approach can be employed to calculate the size distribution and further utilized in modeling, the pellet size distributions are considered using the volume approach while simulating sample beds T1 and T2, as they are closest to the manufacturers' data, as depicted in Fig. 9.
The volumetric approach also allowed us to gauge the average bed porosity using: Based on the water volume, which includes the void volume between the pellets and the total volume occupied by both pellets and water, the average bed porosity was determined to be 45.5% (T1) and 47.1% (T2).These results seem at first close to the 50% porosity for spherical particles used by Beheshti et al. [39].The non-sphericity of our pellets, however, results in a denser bed structure, with the lowest bed porosity occurring for the pellets with the lowest circularity.Table 9 shows good agreement between the bulk bed porosity of the DEM-generated beds and the experimental data.A possible explanation for the deviation is that the DEM packing is not perfectly dense, in the sense that the container is shaken during the experimental construction of the bed.This can be observed from the bed height difference, with the DEM-generated bed for T1 being about 1 cm higher than the real one.

Radial averaged bed porosity
The 3D structure of the solid bed for sample T1 was visualized using X-ray tomography.As depicted in Fig. 10, the irregular shape of the pellets results in a less uniformly structured bed.This is evident in a vertical cross-section of the bed along the y-axis (Fig. 10 (a)), which shows the presence of larger voids distributed across the bed alongside localized denser regions.
The ratio of each pixel belonging to a pellet or void space is collected along the diameter of the bed.An example of the resulting porosity profile for a x-y slice through the bed is shown in Fig. 10 (b).The periodic minima in the first 2 cm from the wall indicate some underlying order in the bed, particularly along the wall sections.This order is often broken up in the middle of the bed, as seen in the x-y slice in Fig. 10 (a), resulting in less pronounced minima in the porosity profile.The structure of a fixed bed reactor is typically characterized by the local radialaveraged bed porosity and bulk bed porosity.
In order to gain confidence in the DEM simulations, these characteristics are compared to the experimental data and an empirical correlation.In the modeling, the local radial bed porosity is calculated as follows: Here, 33 concentric planes of equal thickness are considered in the calculation, as illustrated in Fig. 11 (d).In the above equation, r, A solid , and A total represent the radial distance of the plane from the center, the cross-sectional area of the pellet's in that particular plane, and the total area of that particular plane, including the solid and void areas, respectively.As we can see in Fig. 11 (a), the DEM simulation curve looks like a damping oscillatory function, with the first minimum at 0.5, followed by the first maximum at 1. Near the wall, the radial axial porosity is 1.For the bed T1 in Fig. 11 (a), the DEM-simulated radialaveraged bed porosity profile closely follows the periodicity of the de Klerk correlation [13] and the current experiment, with acceptable deviations in magnitude towards the center.The de Klerk correlation is an improved radial bed porosity and averaged bed porosity prediction model proposed by Arno de Klerk, which considers equal-sized spheres  to address the drawbacks of the existing radial bed porosity models.The deviations of the DEM data with the current experiment increase towards the center of the bed due to the increasing number of irregular/non-spherical pellets accumulated in the center region of the bed, as seen in Fig. 10 (b).These deviations are therefore explained from the fact that the experimental pellets are not spherical, while the DEM simulations are performed on spherical pellets with variable diameters.The mean circularity of the pellets in Fig. 11

Dynamic flow and reduction behaviors in the beds
The gas velocity profiles at 500 sec.are shown in Fig. 12 (a) for all three bed structures in a planar contour plot and exhibit a maximum axial velocity of 5.9 times the superficial inlet velocity (~1 m s − 1 ).To ensure a uniform flow field for the reducing gas, an empty bottom layer has been introduced.The empty bed at the bottom of the pellets helps to minimize the backflow effect.The gas velocity distribution in exemplary slices of T1 DEM , T2 DEM and T1 CT shows significant differences.But in all of these slices, higher magnitudes of velocity are visible in the narrower spaces, and some stagnant zones are scattered across the plane.The gas velocity increases as it leaves the bed, particularly in the top part of the bed slices.The overall gas velocity distribution inside the T1 CT appears to be higher and more uniform than in T1 DEM and T2 DEM .This is also evident from the uniform H 2 concentration propagation inside the T1 CT in Fig. 12 (b).Furthermore, the gas velocity gradually increases with the height of the reactor bed.The H 2 concentration is higher in the lower part than at the top part of the bed due to the flow direction and the reaction progress, as we can see in Fig. 12 (b).The high concentration of H 2 gas in the bottom layer of the reactor results in increased reactant presence inside the pellets in this region and a higher reduction rate.The lower left and the upper right sides of the T2 DEM bed display wall effects on gas velocity, with velocities up to ~6 m s − 1 .The channeling effects in  indicate a lower reduction rate.But this is only a single plane, therefore more quantitative results are given in the following.The gas velocity and hydrogen fields are given at other time steps in the supplementary material.The effects of pellet size and shape on flow dynamics can be summarized as follows: the bed with larger pellets tends to have larger voids, causing uneven flow distribution and potential channeling.In contrast, smaller pellets create smaller voids, resulting in a more even flow distribution throughout the bed.Spherical particles, due to their symmetrical shape, generally promote uniform flow, reducing localized high-velocity regions and enhancing overall flow distribution.Conversely, irregularly shaped particles increase flow resistance and create complex flow paths, leading to higher pressure drops and more stagnant zones.
In Fig. 13, the CFD results for the global conversion degree (F) of the three beds T1 DEM , T2 DEM , and T1 CT are plotted.Since the reduction reactions of iron ore with H 2 are endothermic in nature, the bed temperature will deviate from the initial temperature.All three simulations account for this by incorporating the experimental bed temperature, as shown in the subplot of Fig. 13.The conversion degree plots clearly show that the reduction rate of T1 CT is very close to the one of T1 DEM , and that T2 DEM is the slowest of all.The close agreement observed between the T1 DEM and T1 CT indicates that, despite some discrepancies in the bed morphology between the actual structure and the numerically generated structure, the DEM model is sufficiently reliable.
Table 10 compares the reduction efficiency and rate of reduction for two different bed structures, T1 and T2, at various conversion degrees (50%, 75%, and 95%).The T1 bed structure consistently achieves faster conversion times than T2 across all stages.Specifically, T1 reaches 50% conversion in 506 sec.And 95% conversion in 1522 sec., while T2 takes 680 sec.And 2230 sec., respectively.The rate of reduction is also higher for T1, particularly between 50% and 75% conversion, where it achieves a rate of 0.064 %.sec − 1 compared to 0.043 %.sec − 1 for T2.
Fig. 14 gives, for the same plane as Fig. 12, the local conversion degree for all three beds at 500 and 1000 sec.The conversion degree at 1500 sec.can also be seen in the supplementary material.As discussed in the first part, the plane is not cutting the pellets equally close to their centers, so comparison needs to be done on pellet cross-sections of similar areas.As for the hydrogen flow fields, one cannot draw generalized conclusions, as it is only one exemplary plane for each bed among others.But from Fig. 14, the wall effect seems particularly important for   T2 DEM and T1 CT than T1 DEM , with much further reduced pellets close to the walls.The bottom layer of T1 CT is also clearly more reduced than the pellets on the top and near the central axis.On the contrary, the reduction degree of the pellets in the T1 DEM simulation appears more uniform across the bed.As it could be expected, the overall degree of reduction of the large pellets of T2 DEM is the lowest, featuring pellets with unreacted cores.This indicates an internal diffusion limited process.Fig. 15 shows the averaged local conversion degree for three bed structures at 500, 1000, and 1500 sec., from the wall towards the center (a-c) and from the bed bottom to the top (d-f).These results confirm the former conclusions.For all three beds, the state of reduction is highest at the wall vicinity and globally decreases along the axial direction.There is a visible exception to this axial decreasing trend for T2 DEM , which presents a minimum around 8 mm at all time steps.It can be explained as follow: The bottom pellet layer is well-arranged and 8 mm corresponds more or less to the position of the first pellets' center, while above it, the cross sections meet centers as well as edges.Overall, both the T1 DEM and T1 CT are quite close in magnitude and periodicity.Divergence is higher towards the center of the bed.This happens due to the difference in the radial bed porosity in the central part of the bed (Fig. 11 (a)).We retrieve the overall lower reduction degree of case T2 DEM , but we also see a larger amplitude.This might be due to the lower number of pellets on each section, amplifying the local variability.

Conclusion
This research work utilizes experimental construction of iron ore beds, computational reconstruction via the discrete element method and via CT-scan data, followed by reactive CFD simulations, to evaluate the effect of bed structures on the overall conversion of iron ore pellets for the H 2 -based reduction in fixed bed reactors.The findings highlight that reactor operators should prioritize the selection of pellets with uniform size and spherical shapes, as these configurations enhance gas flow uniformity and narrower reduction rates.In regards of the size distributions of the two pellet samples, T1 and T2, the minor axis of an ellipse fit from image analysis and the volumetric approach based on the water displacement method gave very similar results.The pellet area analysis, however, overestimated the diameter, because the 3D structure of the pellet was largely ignored.
The T1 DEM and T2 DEM beds are computationally reconstructed using the size distributions from the volumetric approach, while T1 CT is a direct reconstruction using the 3D tomographic image.The radial average bed porosity and bulk bed porosity of these DEM-reconstructed beds are compared against experimental and empirical data.The DEMcomputed beds appear loosely packed compared to the real cases.Deviations in bed height and radial average bed porosity can be explained from the limitation of the DEM simulations to use spherical pellets, while small, irregular pellets piled up in the center of the bed in the actual experimental bed construction.Shaking of the vessel further accentuated the real bed packing.This finding suggests that reactor designs should incorporate flexibility in bed packing and porosity adjustments to accommodate different pellet characteristics.
A mesh generator was then used to generate a mesh for each bed with pellet domains refined enough to capture the gas and solid species evolution during the reduction process.Comparison between CFD flow fields, H 2 concentration profiles, local and global reduction degrees has permitted to gain more insight on the influence of the particle size and shape on the reduction process.It could be evidenced that, for beds of similar weight, the bed with a larger average pellet size reduces slower than the bed with a higher number but smaller pellets.This implies that reactor operators may need to dynamically adjust gas flow rates based on real-time monitoring (e.g.pressure sensors, temperature probes, flow meters etc.) of bed structure to optimize efficiency.The findings also suggest that different bed structures may require tailored operating conditions to maximize efficiency.For instance, beds with higher porosity might benefit from adjusted gas flow rates to ensure even distribution and reaction kinetics.The simulation of T2 DEM also showed strong wall effects and channeling, as well as a higher local reduction degree amplitude.Despite the differences in the bed structures between T1 DEM and T1 CT , the global conversions over time are extremely close.There might be several effects counter balancing, for example higher surface areas in T1 CT case but more uniform hydrogen flow across the bed T1 DEM .But the overall good agreement between the cases suggest that the DEM-numerically generated bed is representative enough to model the reduction process of a real fixed bed process when equivalent total weight and pellet size distribution are used.
Based on these findings, it can be also be pointed out that the optimization of a fixed-bed reactor for the direct reduction of iron ore pellets requires a narrower pellet size distribution.A sample containing a wide range of pellet sizes poses practical challenges for reactor optimization.The presence of high variations in pellet sizes may result in incomplete reduction for certain pellets or require increased residence time and gas consumption to achieve complete reduction.This implies that a more uniform pellet size distribution is critical for efficient and practical reactor operation in iron ore pellet reduction processes.Incorporating these insights into reactor operation and design will enable more efficient and effective reduction processes in industrial settings, ultimately enhancing overall reactor performance.
Moving forward, the importance of investigating spatial and temporal temperature variations is recognized, as these factors significantly impact reduction efficiency.Future research will extend current modeling efforts to include other reducing gases, such as syngas and CO, to assess their effectiveness and optimize the reduction process.Additionally, a scaled-up model is being explored to examine how findings from smaller-scale fixed-bed reactors translate to larger industrial applications.Furthermore, the effects of varying gas flow rates will be studied, providing further insights into optimizing reactor design and operation under different industrial conditions.

Fig. 2 .
Fig. 2. (a) Tomographic image capturing of a fixed bed; (b) 3D bed slicing technique to calculate the average radial bed porosity; (c) Axial slices considered for radial porosity calculation at different geometrical degree.

Fig. 3 .
Fig. 3. Illustration of packed bed generation by using DEM simulations at different time steps (a), (b), (c), and (d) generated bed geometry using the input from DEM simulations when all the particles are settled.

Fig. 4 .Fig. 5 .
Fig. 4. Mesh structures of the T1 DEM , T2 DEM and T1 CT along an exemplary plane.The highlighted areas are intended for better visualization.The rightmost figures illustrate the gas and solid-phase mesh structures in real bed.

Fig. 6 .
Fig. 6.Schematic of the iron ore pellet packing model and boundary conditions used for the CFD simulations.

Fig. 7 .Fig. 8 .
Fig. 7. Example image analysis of the samples.(a) Overhead image of the pellets, (b) binary image, (c) outline of the fitted area, (d) ellipse fit corresponding to each pellet sample.

Fig. 9 .
Fig. 9. Size distribution of samples T1 and T2 from image analysis based on the diameter extracted from the minor axis of an ellipse fitted to each pellet (blue), and the diameter extracted from the assumed spherical volume of each pellet (red).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) (c) shows this uneven distribution of the pellets.After the first maximum at 1, the pellet circularity drops, and the lower circularity values indicate the irregular pellets in the center region of the bed.The DEM simulation results for the T2 sample are in line with the correlation data shown in Fig. 11 (b).

Fig. 10 .
Fig. 10.(a) An exemplary slice along the horizontal and vertical axes through the center of the solid bed, and (b) the corresponding averaged porosity profile of the single slice along the bed diameter.

Fig. 11 .
Fig. 11.(a), (b) Comparison of the radial averaged bed porosity distribution for T1 and T2 bed structures against the de Klerk correlation [13], and current experimental data; R, r and d represent the radius of the reactor, the radius of the corresponding concentric plane and the average diameter of the pellets, respectively.The average diameters considered for T1 and T2 are 11.8 and 15.0 mm, respectively; (c) The mean circularity distribution of the real bed; (d) cylindrical planes for the calculation of the radially averaged porosity.

Fig. 13 .Fig. 14 .
Fig. 13.Comparison of the global conversion degree over time for T1 DEM , T2 DEM , and T1 CT respectively.The experimental temporal temperature profile in the inset is taken from Beheshti et al. [39].

Fig. 12 .
Fig. 12.(a) Gas velocity contour at 500 sec. in the T1 DEM , T2 DEM , and T1 CT bed structures.The planer contour plots refer to the y-component of the fluid velocity, where y is parallel to the main flow direction; (b) H 2 gas concentration profiles for the corresponding structures.

Fig. 15 .
Fig. 15.(a), (b), (c) Comparison of the evolution of the radial distribution of the local conversion degree for T1 DEM , T2 DEM , and T1 CT at times 500, 1000 and 1500 sec., respectively; (d), (e), (f) Comparison of the evolution of the axial distribution of the local conversion degree for T1 DEM , T2 DEM , and T1 CT at times 500, 1000 and 1500 sec., respectively.

Table 2
Characteristics of iron ore pellet samples used in DEM simulations.

Table 3
Modeling parameters used in DEM simulations.

Table 4
Coefficients parameters used in collision models.

Table 5
Mesh cell counts.

Table 6
Boundary conditions of the CFD simulation.

Table 7
Surface reaction steps and the kinetic rate parameters of iron ore direct reduction process.

Table 8
Comparison of total pellet volumes from each fitting approach.

Table 9
Comparison of bulk bed porosity of DEM-simulated beds.

Table 10
Comparative reduction efficiency and rate of reduction for different bed structure.