Gianluca Boccardo

Frederic Augier

Yacine Haroun

Daniel Ferre

Daniele L Marchisio

Validation of a novel open-source work-flow for the simulation of packed-bed reactors

Keywords: packed-bed reactor, catalytic reactor, Blender, OpenFOAM, Ergun law, scalar dispersion

The simulation of flow and transport in packed-bed (catalytic and non-catalytic) reactors is of paramount importance in the chemical industry. Different tools have been developed in the last decades and most of them relay on Discrete Element Method (DEM), for generating the particle packing, and Computational Fluid Dynamics (CFD), for simulating fluid flow and scalar dispersion. This work-flow presents the main drawbacks of being computationally expensive, as DEM codes are designed to describe with very high accuracy particle-fluid interactions, that very often are negligible during packing generation, of dealing with non-convex objects, such as trilobes, with cumbersome strategies, and of making use of in-house or commercial codes, that are either difficult to access or costly. In this paper a novel open-source and easily accessible work-flow based on Blender, a rigid-body simulation tool developed for computer graphics applications, and OpenFOAM, a very well-known CFD code, is presented. The approach, which presents the main advantage of being computationally fast, is validated by comparison with experimental data for global bulk porosity, local porosity distribution and pressure drop. To our knowledge this is the very first application of Blender for the simulation of packed-bed reactors.

Introduction

In chemical engineering many processes involve packed-bed (mainly catalytic) reactors and as such the study of momentum, heat and mass transport inside these systems is of great interest. Their investigation is generally performed by using both experimental and modeling techniques. The development of mathematical models, as in other areas of chemical engineering, is particularly interesting because it allows to collect more information and gain more insight, making it easier for scientists and engineers working in this field to investigate the relationship between (catalytic) particle properties (i.e. shape and size distribution) and packing features [START_REF] Lopes | Assessment of cfd euler-euler method for trickle-bed reactor modelling in the catalytic wet oxidation of phenolic wastewaters[END_REF][START_REF] Shahamiri | Modeling catalytic oxidation of lean mixtures of methane-air in a packedbed reactor[END_REF][START_REF] Lopes | Three-dimensional numerical simulation of pressure drop and liquid holdup for high-pressure trickle-bed reactor[END_REF][START_REF] Jafari | Modeling and cfd simulation of flow behavior and dispersivity through randomly packed bed reactors[END_REF]. Moreover, after validation with experiments the working model can be efficiently used to design, scale up and optimize such systems, without resorting to costly and time-consuming experimental campaigns [START_REF] Behnam | A new approach to fixed bed radial heat transfer modeling using velocity fields from computational fluid dynamics simulations[END_REF][START_REF] Dixon | Flow, transport, and reaction interactions in shaped cylindrical particles for steam methane reforming[END_REF][START_REF] Behnam | Comparison of cfd simulations to experiment under methane steam reforming reacting conditions[END_REF][START_REF] Dixon | Experimental validation of high reynolds number cfd simulations of heat transfer in a pilot-scale fixed bed tube[END_REF][START_REF] Romkes | Cfd modelling and experimental validation of particle-to-fluid mass and heat transfer in a packed bed at very low channel to particle diameter ratio[END_REF]. An additional advantage of developing mathematical models in this field, is that the mathematical model can be used to identify numerous parameters that characterize the fixed-bed, or packing, such as porosity and permeability, that might, under some operating conditions, be hard to determine experimentally [START_REF] Augier | Numerical simulations of transfer and transport properties inside packed beds of spherical particles[END_REF][START_REF] Augier | Numerical approach to predict wetting and catalyst efficiencies inside trickle bed reactors[END_REF].

A mathematical model for packed-bed reactors is mainly constituted by two modules: one that generates the geometry of the packing and one that simulates fluid flow. Usually the second task is performed by running Computational Fluid Dynamics (CFD) simulations [START_REF] Atmakidis | Cfd-based analysis of the wall effect on the pressure drop in packed beds with moderate tube/particle diameter ratios in the laminar flow regime[END_REF][START_REF] Mousazadeh | Eulerian simulation of heat transfer in a trickle bed reactor with constant wall temperature[END_REF][START_REF] Shi | Cfd analysis of liquid phase flow in a rotating packed bed reactor[END_REF][START_REF] Dixon | Cfd method to couple three-dimensional transport and reaction inside catalyst particles to the fixed bed flow field[END_REF][START_REF] Dixon | Fixed bed catalytic reactor modelling-the radial heat transfer problem[END_REF], whereas for the first task two strategies can be used. The first lies on using actual experimental data, be it via Scanning Electron Microscope (SEM) imaging of the catalytic particles [START_REF] Boccardo | Microscale simulation of particle deposition in porous media[END_REF], micro-computer tomography scans of a portion of the packed-bed loading [START_REF] Mutailipu | Pore network simulation of two phase flow based on x-ray micro computed tomography images[END_REF][START_REF] Akolkar | Tomography-based characterization and optimization of fluid flow through porous media[END_REF][START_REF] Schmidt | Estimation of chromatographic columns performances using computer tomography and cfd simulations[END_REF], or other such methods. While these procedures have a merit in the precision with which the description of the medium at the scale of the voids between the particles is obtained, they suffer from a great difficulty in the post-processing of these scans, in order to extract a suitable mesh for the CFD code. The second possibility, is to rely on an in-silico algorithmic reconstruction of a geometrical model, which faithfully represents the real system in all of its features: the first advantage is the extremely lower cost of such an approach with respect to using the more sophisticated experimental techniques mentioned, coupled with the easiness of generating and testing a very high number of loading realizations, with varying particle shapes and particle size distributions. Obviously, great care must be given as to ensure that the reconstructed model is accurate both as a purely geometrical description and, perhaps more importantly, in showing the same fluid dynamic behavior of the real system under investigation.

The in-silico algorithmic reconstruction is often performed by using the Discrete Element Methods (DEM) [START_REF] Horgue | Experimental and numerical study of twophase flows in arrays of cylinders[END_REF][START_REF] Tabib | A 3d cfd-dem methodology for simulating industrial scale packed bed chemical looping combustion reactors[END_REF][START_REF] Vollmari | Pressure drop investigations in packings of arbitrary shaped particles[END_REF][START_REF] Mahmoudi | Application of xdem as a novel approach to predict drying of a packed bed[END_REF][START_REF] Tsory | Thermal dem-cfd modeling and simulation of heat transfer through packed bed[END_REF][START_REF] Bai | A coupled dem and cfd simulation of flow field and pressure drop in fixed bed reactor with randomly packed catalyst particles[END_REF] that present the main disadvantages of being computationally expensive, because of the accuracy with which fluid-particle interactions, negligible during packing generation, are described and being very often limited to the treatment of simple convex particles (i.e. spheres), whilst more realistic particle shapes (i.e. cylinders and trilobes) can be treated only with cumbersome strategies. The objective of this work is therefore to introduce, test and validate an innovative approach for packing generation, alternative to DEM, still based on the classical ballistic theory of rigid-body simulations, but implemented in the very fast and computationally efficient computer graphics code Blender [START_REF] Van Gumster | Blender For Dummies[END_REF]. An advantage of Blender is the capability of dealing with particles of non-convex shapes, such as trilobes, very much used in catalytic reactors. The geometry in this way generated is than used to build a computational mesh and run CFD simulations in the OpenFOAM environment [START_REF] Opencfd | The Open Source CFD Toolbox, User Guide[END_REF]. It is indeed interesting to observe that the developed work-flow, from geometry and mesh generation, to CFD simulation and post-processing totally relies on free open-source codes.

Since to our knowledge, this is the first application of the proposed open-source workflow to the simulation of packed-bed reactors, a full validation is required. Therefore the capability of the tools to build realistic packings for three different catalyst particle shapes, namely spheres, cylinders and trilobes, is assessed, by comparing our predictions with experimentally measured packing bulk porosities, as well as local radial porosity profiles in cylindrical vessels. Subsequently, for each catalyst shape, the influence of some features of the particle size distribution on the resulting packing bulk porosity is analyzed. Also, a great deal of attention is paid to the meshing process. Meshes with different levels of refinement are tested and results concerning the evaluation of the critical cell size necessary to obtain both a good discretization of the particle surface and a precise description of the contact points between them, are discussed. CFD simulations are eventually performed on these systems by solving the continuity and Navier-Stokes equations and predictions regarding pressure drops in the bed are validated by comparison with empirical correlations. The final results prove that the proposed open-source work-flow is capable of simulating packed-bed reactors with very high accuracy and low computational costs.

Governing equations and theoretical background

Blender solves the Newton's equations of motion for a system constituted by N particles or bodies of arbitrary shape. It does that by availing itself of the Bullet Physics Library (BPL), which is a large collection of codes used to manage the dynamics of rigid-bodies and, most importantly, to detect and calculate the outcome of the collisions between these bodies. This library provides for a number of iterative methods combining accuracy, speed and robustness, enabling for the simulations of a very large number of elements, as it will be shown further on; moreover, a clear advantage of using this code with respect to many other algorithms used to recreate granular media models, lies in the possibility to manage any particle shape, even complex non-convex ones. An additional advantage lies with the fact that Blender does not consider the flow of the fluid inbetween particles; being generally this fluid air, its presence has very little influence during the packing process. This makes the use of Blender more interesting, from the computational point of view, with respect to other similar codes, that instead do account for these effects (i.e. DEM). Details on the BPL can be found in the literature [START_REF] Coumans | Bullet physics library[END_REF].

Once the geometry is generated, it is discretized and meshed within the OpenFOAM environment, which is an open-source finite volume CFD code. The equations solved concern with the momentum and mass balance at the pore-scale (i.e. the scale of the voids left between the particles). For a fluid of constant density and viscosity, these are the well known Navier-Stokes and continuity equations; for a quick summary of these pore-scale equations, and their relationship with macro-scale equations, we remand to our previous work [START_REF] Icardi | Pore-scale simulation of fluid flow and solute dispersion in three-dimensional porous media[END_REF]. Here we limit the discussion to the macro-scale equations used for validation in this work.

From the macro-scale point of view, the pressure drop per unit length, ∆P/L, in packed-bed reactors is calculated by using Ergun's law:

∆P L ρD g G 2 0 ε 3 (1 -ε) = 150 1 -ε (D g G 0 )/µ + 1.75, (1)
where G 0 = ρq is the mass flux, ρ is the fluid density, q is the superficial velocity, ε is the bulk porosity, µ is the fluid dynamic viscosity and D g is the size of the particles constituting the packing. D g is generally taken equal to the actual diameter in the case of spheres and in the case of circular cylinders with equal diameter and height, whereas different definitions of equivalent diameters are used for other objects. For example, for trilobes it is obtained via the empirical correlation: D g = CD t , where D t is the circumscribing diameter [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Twophase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF]. The empirical factor C depends on the trilobe size, and was found to be equal to 0.8 for the trilobes investigated in this work.

Test cases and numerical details

In the following section the different test cases explored, operating conditions employed and numerical detailed involved in this work are described, starting from the generation of the geometric models, down to their meshing and concluding with the fluid flow simulations.

Obtaining the geometric model. As mentioned, the first part of this work deals with the development of a robust tool for effectively simulating the stacking of particles in a packed bed, with particular attention to catalytic particles and catalytic reactors. In order to do that the software Blender [START_REF] Van Gumster | Blender For Dummies[END_REF] is used. The rationale of this choice, besides the arguments already mentioned, lies also in its extensive scripting functionality. A tool was then developed (in the language Python) and plugged into the main rigidbody simulation code, with the purpose of quickly setting up different simulation cases, with several particle shapes, particle size distributions, total numbers, and containers in which the particles settle. The shapes considered are spherical, cylindrical and trilobe extrudates, while the containers tested were cubes and cylinders. An example of three loadings of spheres, cylinders, and trilobes (in a cylindrical container), generated with Blender can be found in Fig. 1.

Since the purpose is to mimic the loading of the catalyst beads into the reactor in a realistic manner, the particles are initially placed on top of the container, with gravity causing their deposition into it, with the end state of the simulation provided by a stable solution of the balance of the forces acting on them, namely gravity and the interaction forces. The choice of the initial placement of the grains in the setup of the rigid-body simulation is a step worth discussing. Three different strategies are tested. In short, the first way is the simplest to implement but the less realistic one, with every particle falling in the container through a single straight line, while the other two represent the actual process quite well, with the random "grid" placement giving the best results. The difference in the latter lies in the additional step of allowing for a certain distance between the grids, resulting in more time for each set of particles to settle and thus in a more stable simulation. For non-spherical items, such as cylindrical beads and trilobes, a randomization of their orientation is also added during their placement in the grid.

Another very important point in the setup concerns the geometric description of the particles. All particle shapes are represented in Blender with "watertight" (closed) external meshes characterized by different refinement levels. For example, a sphere with refinement level equal to one is represented as an icosahedron, with each further refinement level splitting each triangular face into four triangles, resulting in a smoother surface.

Particle-particle interactions are another very important phenomena involved during packing generation. In our work a hard-sphere potential was employed, implying no interpenetration, and transforming interactions into instantaneous collisions. These are computed in Blender with different strategies. In the case of spherical particles collision is simply defined through the centers of the colliding spheres and their radii. In the case of more complex shapes, the computation of the nearest-approximating convex hull is needed, together with the surface mesh, leading to higher computational costs. Also the coefficient of restitution (identifying the ratio between the post-collision and the precollision kinetic energy) needs to be specified. Besides collisions also enduring contact (i.e. sliding) between particles is considered; this latter is governed by the friction factor, another parameter that need to be specified.

Both of these parameters (i.e. restitution coefficient and friction factor) can be modified, greatly influencing the outcome of the physical simulation. For example values of the restitution coefficient closer to unity generally lead to denser packings, as particles are left, collision after collision, with larger fractions of their pre-collision kinetic energy, resulting in more persistent vibrations and further adjustment of particle relative positions. The same effect is obtained by reducing the friction factor. The final values employed are often the result of parameter identification performed with the objective of improving the agreement with experimentally measured packing bulk porosities [START_REF] Augier | Numerical simulations of transfer and transport properties inside packed beds of spherical particles[END_REF].

In the results section the optimal choice for all these parameters (i.e. refinement level, coefficient of restitution and friction factor) will be presented and justified. This choice was based on the investigation of three different packings. The first one was constituted by spheres with average diameter equal to 1.99 mm, the second by cylindrical beads with both average height and diameter equal to 1.31 mm and the third by trilobes, with circumscribing average diameter of 1.8 mm (and equivalent diameter equal to 1.44 mm), and length equal to 4.77 mm. Both perfectly mono-disperse and poly-disperse populations of spheres, cylinders and trilobes were considered. The average diameter reported above refers to the actual diameter of the objects in the mono-disperse case, whereas it corresponds to the mean diameter of a Gaussian distribution with variance σ 2 in the poly-disperse case; the standard deviation σ (square root of the variance) was varied roughly in between 0 and 0.9 mm. Packings were generated in triplicates in cylindrical containers of equal diameter and height, with these values being equal to 3.15 cm, 3.57 cm and 4.3 cm respectively for spheres, cylindrical beads and trilobes packings.

Mesh Generation. A mesh is thus created for these three packings, to be used in the finite-volume CFD code, for the fluid flow and scalar dispersion simulations. For these simulations cubic samples (with sides equal to 1.25 cm) extracted from the larger cylindrical containers described above are employed. The size of these samples is accurately chosen in order to have a representative elementary volume (REV) as detailed later on. The mesh utility snappyHexMesh, included in the open-source package OpenFOAM [START_REF] Opencfd | The Open Source CFD Toolbox, User Guide[END_REF], is used at this end. The first step in the meshing process is the creation of a structured, cartesian grid in the entirety of the bounding volume of the system considered. Next, the volume pertaining to the geometric model is subtracted, leaving a cartesian grid in the fluid portion approximating the particles surface with a stair-step description. Finally, the mesh is further modified by moving each boundary vertex, relocating them closer to the original model surface, resulting in a body fitted mesh. The mesh generation process is critical to both obtaining grid independent results during the fluid flow simulations (i.e. resolution of the boundary layer around particle surface), and for ensuring that the model constitutes a faithful representation of the original geometry (i.e. definition of the contact points), as described in the literature [START_REF] Dixon | Systematic mesh development for 3d cfd simulation of fixed beds: Contact points study[END_REF].

Examples of the final models after the meshing operation, for the three particle shapes considered (i.e. spheres, cylinders and trilobes) are reported in Fig. 2, while Fig. 3 shows a close-up view of the actual mesh for one of them. Each of the samples is comprised of approximately 300 grains, with the size of the mesh being respectively equal to around 8.5, 11 and 22 million cells. Table 1 explains the meshing strategy employed in each case, showing the two different cell dimensions for the "near-grain" zone (four cells thick for all cases) and the volume bulk.

Fluid flow simulations. Using the meshes on the cubic samples obtained as described in the previous part, fluid dynamic simulations are performed by using the CFD code OpenFOAM 2.3.0 (some simulations were also run with the code Ansys Fluent 15 for comparison purposes but are omitted for brevity). The fluid velocity field is determined by solving the continuity and Navier-Stokes equations, after imposing a given pressure drop. A no-slip boundary condition is applied on the grain surface. An inlet zone is set on one side of the geometry and an outlet zone at the opposite face, with the fluid entering the domain due to an imposed pressure drop through the medium. No additional forces (like gravity), are considered. On the remaining sides a condition of symmetry is set, implying no fluid flow across these surfaces, as it is shown in Fig. 4. The Newtonian incompressible fluid considered here is water, with density ρ = 997.78 kg m -3 , dynamic viscosity µ = 9.77 × 10 -4 kg m -1 s -1 . An overview of the operating conditions of each case for each grain shape can be found in Tab. 2 and as it is seen the Reynolds number ranges from about 10 -3 to 10 2 . The system is solved for in laminar conditions and under constant temperature (i.e. the energy equation is not solved), in steady-state conditions. Flow field data extracted from the simulations are compared to the predictions from the well-known Ergun law, reported in Eq. [START_REF] Lopes | Assessment of cfd euler-euler method for trickle-bed reactor modelling in the catalytic wet oxidation of phenolic wastewaters[END_REF]. It is worth mentioning that the mean particle diameter appearing in the Ergun equation, D g , is taken equal to the actual diameter in the case of spheres and cylindrical beads, while for the trilobes it is the calculated equivalent diameter: D g = CD t , where D t is the circumscribing diameter and As already introduced when discussing the necessity of properly resolving the momentum boundary layer around each particle, grid independence of the CFD results also needed to be ensured, and the parameter used to assess this is the equivalent particle diameter D * g . Following a methodology used and described in greater detail in a previous work [START_REF] Boccardo | Microscale simulation of particle deposition in porous media[END_REF], each different mesh refinement level is characterized by an effective diameter D * g , which is identified by fitting the pressure drop results from steady-state flow simulations at different fluid velocities with the Ergun law. Variations of D * g versus total number of mesh cells (or cell dimension) can then be analyzed and the mesh ensuring grid independent results identified.

Spheres

Results and discussion

As mentioned many times our aim is to develop a comprehensive work-flow for the simulation of packed-bed reactors with open-source simulation tools. In order to do that, the realisticness of the geometrical models obtained and of their fluid dynamic features have to be assessed. In this section our results are presented, starting from how a representative elementary volume of the considered media is obtained to how the geometrical features of the models (bulk and radial porosity) compare with the experimental data.

Moreover, the effect of the distribution variance on the resulting packing bulk porosity is also analyzed. Next, we perform simulations of fluid flow, under different operating conditions, in the different types of packing prepared.

Geometrical models. As anticipated, the first test has the purpose of identifying a suitable "size" for the simulation (in terms of the total number of particles), in order to quantify the dimension of the smallest REV for the system under investigation. This analysis was performed on different types of packings, but for the sake of brevity only results concerning one case will be presented here. A cubic container, 2.3 cm in size, is filled with 1772 spherical particles, distributed along a Gaussian distribution with an average diameter of 1.99 cm and a standard deviation σ = 0.29 mm. Cubic samples of different dimensions are then taken inside the system, all centered in the center of gravity of the domain, and the corresponding porosity is computed. The results are shown in Fig. 5 for three different refinement levels. At smaller dimensions, very high porosities are reported: this is expected as, for sample sizes of the order of the particle diameter, a disproportionately high fraction of a void or a particle could end up in the sample. At 5 mm and for bigger sizes, the reported porosities reach a stationary value corresponding to about 0.355, showing another instability arising only with sample sizes of the order of the container size, where wall effects come into play [START_REF] Zou | The packing of spheres in a cylindrical container: the thickness effect[END_REF]. Finding this behavior for this kind of analysis is very satisfactory and equivalent results can be found in other examples in the literature [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. Moreover, this offers an indication for the choice of a refinement level for the particles description: as it can be seen, there are not substantial differences between the three cases except for small sample sizes. What is indicated as refinement level 3 is the most suitable for these applications, as it combines a good representation of the actual shape with an acceptable computational cost during the packing simulation.

Cylindrical containers are also tested, as they resemble the actual reactor shape and as lab-scale experiments are usually carried out in these test containers. A poly-disperse Gaussian distribution of 4000 particles, with an average diameter equal to 1.99 mm, and a standard deviation of 0.29 mm, is loaded with the random grid technique in a series of different containers, as specified in Tab. 3, and the final bulk porosity is then evaluated. This was calculated by using a dimension of the global sample slightly lower than that of the container itself, to avoid the effect of the presence of the walls, which would affect the final bulk porosity. These values were also compared with those calculated as the average porosities of a series of 50 cylindrical samples, randomly placed in the container (with diameter and height both equal to 5 mm), resulting in very similar values. As expected, since the same particle distribution is used in all cases, the final bulk porosity values are very close to each other and oscillate around ε = 0.360. This result is also in line with data (pertaining to the same particle size) coming from a large experimental campaign dealing with the study of flow conditions in trickle bed reactors [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Twophase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF]. These results were obtained by using 0 and 0.5 for the coefficient of restitution and the friction factor, respectively. As mentioned in the previous section, the final packing characteristics (including bulk porosity) are affected by the values of the restitution coefficient and friction factor, especially in the case of trilobes. As an example Tab. 4 reports the predicted bulk porosity obtained under different settings for these two parameters (for trilobes). The data reported in this table show how strongly these two physical parameters influence the mixing effect mentioned earlier, and hence the resulting bed porosity. Notwithstanding the difficulties in treating the dynamics of rigid-body interactions in the case of such a complex particle shape, it is worth noticing that, when the right set of interaction parameters is used (i.e. cases T2 and T3) a remarkable accordance to experimental data is found [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Twophase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF].

Container

An interesting comparison, that allows to investigate the performance of the packing algorithm is the prediction for the radial porosity profiles, that highlights the well-known wall effects. These were explored for different test cases but for the sake of brevity only a selection of them will be presented. For example, the predictions concerning the packing labeled as P3 (see Tab. 3) are shown in Fig. 6. In this case a different sampling strategy is adopted: 200 annular samples, with height equal to 2 cm, are taken. Each of these was placed at the center of the container, and each has a fixed volume but an increasing diameter, effectively calculating the local porosity value along the radial direction from the center of the packing. In the same figure, these results are compared with two other analyses of the same type: one conducted on an earlier study of an algorithmically generated packing [START_REF] Augier | Numerical simulations of transfer and transport properties inside packed beds of spherical particles[END_REF], and the other one being actual experimental data [START_REF] Giese | Measured and modeled superficial flow profiles in packed beds with liquid flow[END_REF]. The accordance between the results obtained in this work and the experimental data is satisfactory, as it confirms that the physical simulation replicates the characteristics of a real packing consistently, especially in terms of amplitude and period of oscillations.

The same analysis is performed on a packing of uniformly distributed cylinders (with the same geometrical features of the experiment performed in the already mentioned study [START_REF] Giese | Measured and modeled superficial flow profiles in packed beds with liquid flow[END_REF]) and the results are shown in Fig. 7, where they are compared with experimental data. Also in this case good agreement was observed. It is particularly interesting to observe that in the case of cylindrical particles, the experimental oscillations are characterized by spiky maxima and this feature is perfectly reproduced by the packing algorithm.

As last example, the tests conducted with trilobe packings in cylindrical containers are discussed. In this case a Gaussian distribution in length of 3000 trilobes with an average length of 4.77 mm and a circumscribing diameter of D t = 1.8 mm (and equivalent diameter D g = 1.44 mm) and standard deviation equal to 1.8 mm is considered, in a container with diameter and height respectively equal to 3.5 cm and 3 cm. Model predictions are reported for this case in Fig. 8.

Lastly, for each catalytic particle shape, a number of packing realizations (of 3500 grains each) are created, having identical mean particle diameters but differing in the variance of the (Gaussian) size distribution. The features of each of the cases considered are found in Tab. 5, along with the resulting bulk porosity for each packing. It has to be noted that again, as with the analyses reported in Tab. 3, the size of the sample from which the bulk porosity value was obtained is smaller (in the order of two equivalent particle diameters less) than the full container used in the rigid-body simulation, in order to remove the effect of the near-wall region on porosity calculations. Although no experimental data is available for direct comparison, as expected, there is a slight downward trend in porosity with the increase in the particle distribution non-uniformity both in the case of spherical and cylindrical particles, consistently with what found by other authors [START_REF] Farr | Random close packing fractions of lognormal distributions of hard spheres[END_REF]. In the case of trilobes, increases in the distribution variance seem to have little, if any, effect on the packing porosity. This could be due to either an effect specific of the particular shape considered, or to the much different length/diameter ratio of the trilobe particle with respect to the cylindrical ones, causing different long range structures in the packing bringing about this unexpected behavior. This finding calls for deeper investigation, both widening the range of non-uniformity of the distributions and the aspect ratio of the particle shapes considered.

Fluid flow. Before discussing the results concerning this part of the work, it is useful to remind that in all the simulations presented in this work, while building the mesh for simulating the fluid flow, the cell layer next to the surface of the particles underwent an additional refinement step. As mentioned this is done mainly for two reasons. Firstly it is necessary to account for the sharper momentum boundary layer near the walls and thus adequately describing it. Secondly, a certain level of refinement near the surface of the particles is necessary to give a satisfactory definition of the contact points, which can be wrongly approximated as larger contact solid volumes if too coarse a mesh is employed. 13 In order to quantify this effect a 7 mm sized cubic sample was extracted from a larger packing, constituted of a Gaussian distribution of spheres of average diameter equal to 1.9 mm. The porosity of the actual model, result of the Blender simulation, is 0.344. Then, several meshes with different cell dimensions are created, and the ratio between the porosity of the meshed sample and the original geometry, ε m /ε g , is calculated. The results are presented in Fig. 9. As it can be seen, very coarse meshes can result in a very low ε m /ε g ratio, down to ε m /ε g ≈ 0.9, while a grid with average mesh cell size of 50 µm (corresponding to 800,000 cells in this sample geometry) leads to ε m /ε g ≈ 0.99, which corresponds to a faithful representation of the original model. The acceptable cell dimension obviously depends on the average particle diameter: this result shows that an analysis of this type is a necessary step during the mesh generation process. mean diameter in the case of spheres, diameter and average length in the case of cylindrical beads and trilobes. Along, the standard deviation in the distribution (of mean diameter for the spheres and mean length for beads and trilobes) is also reported.

refinement levels. As an example this analysis is reported here for the previous case of a packing of Gaussian distributed spheres. By following the procedure already introduced in our previous work [START_REF] Boccardo | Microscale simulation of particle deposition in porous media[END_REF], instead of plotting the pressure drops at different superficial velocities, data are compactly presented in terms of one single equivalent particle diameter, D * g . As mentioned this is the particle diameter appearing in the Ergun law reported in Eq. ( 1) that results in the pressure drops obtained at different superficial velocities of the fluid flow simulation. Results are reported for this test case in Fig. 10. As it is seen by increasing the number of mesh cells, their size decreases and the corresponding equivalent particle size, D * g , decreases. This happens as by refining the mesh in the momentum boundary layer, the pressure drops are more accurately calculated and result in larger values, corresponding to smaller equivalent particle sizes. As seen, in this case, a cell dimension of 50 µm will result in a satisfactory discretization of the momentum boundary layer as further refinements do not significantly change pressure drops and equivalent particle size. This result is used in the decision of the appropriate cell dimension for the fluid flow simulations for all packings, after having taken into due consideration the different dimensions of the grains and as such the need for a different scale of discretization of the properties boundary layers. As already seen, Tab. 1 reports the final meshing strategy.

Let us now discuss the actual results concerning the fluid flow. Some examples of typical flow fields in the packing of spheres, cylindrical beads and trilobes (for one su- perficial velocity) are shown in Fig. 11. As it is seen, due to the irregular disposition of the particles in the packing, very inhomogeneous flow fields are obtained, with sudden deviations of the fluid and very pronounced velocity profiles within the pores. For all the investigated test cases the observed pressure profiles across the simulated domains was linear. This is a good indicator of the validity of the simulated samples chosen and further confirms that can be geometry can be considered a representative volume.

The aggregate results for all the superficial velocities considered can be found in Fig. 12 for packings of spheres, cylinders and trilobes. As it is shown, the model predictions (symbols) are compared with the predictions of Ergun's law (continuous line), resulting in very good agreement. It is worth mentioning here that the predictions of Ergun's law, reported in Eq. ( 1) and calculated with the particle size, D g , evaluated with the strategy previously described, are considered in this context as experimental data. This is because the predictions of Ergun's law have been extensively (for spheres, cylinders and trilobes) validated on a very large number of experiments. It is therefore legitimate to consider the agreement between the predictions from the model and from Ergun's law as an equivalent validation with experiments, showing once again the realisticness of the generated packings.

Conclusions

In this work a computational tool developed in computer graphics, Blender , integrated with the BPL, is used to generate realistic packings of (catalyst) particles of different shapes. The main advantages of this approach (versus other alternatives) stand in the possibility of simulating packings constituted by particles with complex shapes (e.g. non-convex objects such as trilobes), in the very low computational costs associated, in its extensive scripting functionality and in its open-source nature.

Results show that attention has to be paid to the strategy with which particles are inserted into the container, as well as to the meshing procedure, carried out here with snappyHexMesh. The mesh has to be fine enough to describe well the geometrical details of the pores within the packing, as well as the momentum boundary layer around each particle. The validation with experimental data demonstrates that the generated packings realistically describe the behaviour of catalytic fixed bed reactors.

In particular, when working with poly-disperse spherical particles, the well known radial porosity profiles are obtained, with an accuracy superior to other similar tools. Interesting results are also obtained regarding the influence of the particle poly-dispersity on the resulting packing bulk porosity: while for spherical and cylindrical particles the expected inversely proportional relationship is found, this is not true in the case of trilobes. This calls for a deeper investigation, both with regards to the particle shapes considered and to the particle aspect ratio. This last feature especially could be the cause of unexpected packing structures affecting the final bulk porosity value of the packing.

Moreover, fluid flow simulations, carried out with OpenFOAM 2.3.0 show very good agreement with the predictions of the Ergun law. In particular, when working with polydisperse trilobes very good results are also obtained, especially considering that these complex non-convex objects are very difficult to treat. Also in this case comparison between simulated pressure drops and those predicted with the Ergun law (by using a particle size in turn obtained from corresponding experiments) shows an excellent agreement.

Overall these results show that the open-source work-flow proposed in this work can be used to investigate fluid flow in fixed bed reactors. In the next steps of our work chemical reactions, heat transfer and multi-phase flows will also be considered.

Figure 1 :

 1 Figure 1: Packing of spheres, cylindrical beads, and trilobes in cylindrical containers.

Figure 2 :

 2 Figure 2: Rendering of the final meshed models, respectively for the packing of spheres, cylindrical beads, and trilobes.

Figure 3 :

 3 Figure 3: View of the final mesh, packing of cylindrical beads.

Figure 4 :

 4 Figure 4: Snapshot of the computational domain, with the boundary conditions for the fluid flow in evidence: simmetry boundaries in ochre and the outlet in teal, solid grain walls in grey (inlet zone not pictured). Fluid main direction is indicated by the arrows, red for the inlet, teal for the outlet.

Figure 5 :

 5 Figure 5: Porosity for increasingly large cubic samples for three different spheres refinement levels (respectively equal to two, three and four surface subdivision levels for circles, squares and diamonds).

Figure 6 :

 6 Figure 6: Comparison of local porosity radial profiles: current packing model (white squares), reference packing model [10] (dot-filled circles) and experimental data [35] (black filled diamonds) for a packing of spheres. Distance from wall expressed in number of equivalent diameters Dg.

εFigure 7 :

 7 Figure 7: Comparison of local porosity radial profiles: current work packing model (white squares) and experimental data [35] (black filled diamonds), for a packing of cylindrical beads. Distance from wall expressed in number of equivalent diameters (or equivalently, bead length) Dg.

Figure 8 :

 8 Figure 8: Comparison of local porosity radial profiles: current work packing model (white squares) for a packing of trilobes. Distance from wall expressed in number of equivalent diameters (or equivalently, bead length) Dg.

Figure 9 :

 9 Figure 9: Porosity ratio and number of cells versus cell dimension for a spheres packing.

Figure 10 :

 10 Figure 10: Effective grain size, D * g , with varying number of mesh cells, and corresponding cell dimension.

Figure 11 :

 11 Figure 11: Streamlines and contour plot of fluid velocity (in m s -1 ) in a planar cut for the packing of (from top to bottom) spheres (Re = 1.14), cylindrical beads (Re = 0.28) and trilobes (Re = 0.60). The contour plots refer to the x-component of the fluid velocity, where x is parallel to the main flow direction. Inlet and outlet are on the left and the right side, respectively.

Figure 12 :

 12 Figure 12: Comparison of model results (symbols) with the Ergun's law (continuous line) for packings of (from top to bottom) spheres, cylinders and trilobes.

Table 1 :

 1 Meshing strategy for the three cases considered. Thickness of the "near-grain" layer equal to four refined cells for all cases.

Table 2 :

 2 Operating conditions for each case, divided by geometry considered (particle shape). Cases V1 through V11 differentiate by imposed pressure drop (in Pa m -1 in the table, for each case). In the last row, the range of Re explored for each case, corresponding to the range of pressure drops in the rows above, is reported for each geometry.

Table 3 :

 3 Bulk porosity for spherical particle loadings in four different containers, with the corresponding sample dimensions.

Table 4 :

 4 Restitution coefficient and friction factor used in three different packing simulations of trilobes, with the resulting total and sampled porosity.

		Restitution Friction	ε
	Cases	coefficient	factor	(-)
	T1	0.00	0.500	0.415
	T2	0.70	0.325	0.390
	T3	0.85	0.325	0.378

Table 5 :

 5 The correct numerical resolution of the momentum boundary around each particle in the packing can be assessed by analyzing the predicted pressure drops at different Results of bulk porosity relative to a number of differently distributed packings of spheres, cylindrical beads, and trilobes. Characteristic lengths for each shape are reported:

	Shape	Case	Characteristic	Standard deviation, mm	ε, -
			lengths, mm
		SU1		-	0.3775
		SN1		0.15	0.3720
	Spheres	SN2	D g =1.99	0.22	0.3723
		SN3		0.29	0.3678
		SN4		0.36	0.3652
		CU1		-	0.3444
		CN1	D g =1.31	0.1	0.3407
	Cyl. beads	CN2	L=1.31	0.15	0.3413
		CN3		0.19	0.3410
		CN4		0.24	0.3380
		TU1		-	0.3967
		TN1	D g =1.8	0.36	0.4029
	Trilobes	TN2	L=4.77	0.527	0.4029
		TN3		0.695	0.4037
		TN4		0.86	0.4047

Corresponding author: Daniele Marchisio; Dipartimento Scienza Applicata e Tecnologia, Politecnico di Torino, C.so Duca degli Abruzzi

[START_REF] Mahmoudi | Application of xdem as a novel approach to predict drying of a packed bed[END_REF] 10129 Torino; e-mail: daniele.marchisio@polito.it; tel.: 0039-011-

Acknowledgements

The authors want to express their gratitude to Luigi del Plato for its important contribution during his internship at IFP Energies Nouvelles.