A study of geometrical structure of packed beds using flow simulation with the immersed boundary method

In this work, a new method is proposed for calculation of one of the essentials geometrical parameters characterising a packed bed - tortuosity. The geometry of random packings of cylinders and rings (with various packing densities) is numerically generated with Discrete Element Method approach. For flow calculations in a highly complex geometry of a random packed bed, immersed boundary method is employed (with direct forcing variant). The obtained flow velocity field allows for calculation of tortuosity using two methods: by passive advection of massless particles and by comparison of the mean flow velocity magnitude with the mean streamwise velocity component. It is demonstrated that for the packings of low density these two methods are in a good agreement but may diverge for dense packings. Frequency distribution of tortuosity values for individual flow paths is also examined. The results show more frequent occurrence of low-tortuosity paths for coarse packings.


Introduction
Random packed beds -fixed packings of randomly oriented particles -are widely employed in many areas, especially in chemical and process engineering as a porous medium of a large specific surface enhancing the mass transfer or reaction rate between reacting substances (e.g. in trickle bed reactors or absorbers). The flow characteristics in such beds are thus of a significant practical importance but also are interesting from purely scientific point of view. According to the well-known Darcy formula average flux q of a fluid in the low-Reynolds number regime is proportional to the imposed pressure gradient: where µ -dynamic fluid viscosity and k -permeability of the porous medium directly related to its geometrical structure. In one of the first formulas for k , introduced by Kozeny, permeability depended only on the porosity ε of the medium (or void fraction in the context of packed beds -ratio of volume of empty space between the particles to the total volume of the bed) and its specific surface S (pore surface area per unit volume of the bed) but later was modified by Carman to include another geometrical parameter -tortuosity τ [1,2] where c -structural parameter depending on the cross-section of the pore channels. Tortuosity can be defined in many different ways generally expressing the degree of deformation of flow paths within the medium. For an individual path we define its tortuosity as [2]: where path L is the actual length of the path and L -the thickness (linear size) of the porous medium.
When there are no obstacles in the medium, the path length coincides with L giving 1 τ = but, as in packed beds only the void space can be penetrated by the flow, the paths follow a very complex pattern and 1 τ > (for some structures even 1 τ ≫ ). Tortuosity as a parameter of a medium, not just of an individual path, takes into account an average path length representative to the whole volume of the medium.
In the literature, there are many attempts to find a relation between porosity (void fraction) ε and tortuosity (see a short review in [1]). Such a relation cannot be, of course, universal as it is possible to provide examples of porous media with the same porosity and arbitrary large tortuosities. However, for structures of particular type -e.g. media consisting of randomly distributed grains of a given regular shape -tortuosity can be shown to be a function of porosity through numerical simulations [1,2] or theoretical considerations [3].
One of the approaches to calculation of tortuosity is based on purely geometrical considerations as determination of shortest path between different points of the medium [4,5]; a little different approach is used in [6]. These methods have relatively low computational cost but their main disadvantage is that they do not take into account the physical flow structure in the medium. In principle, it is possible that pore channels comprising the shortest paths are characterised by large flow resistance and the major part of the fluid is transported through totally different paths. In such a case, geometricallybased estimations of tortuosity may be significantly distorted. It seems that physically well-grounded calculations of tortuosity for general packings have to include the information about the real flow.
However, simulation of a fluid flow in such a complex geometrical structure as a random packed bed is definitely not an easy task. Although flow analysis in porous media can be efficiently performed at the level of macroscopic modelling and averaged flow equations [7,8], calculation of tortuosity requires modelling at the microscopic level (of individual bed particles). Application of classical finite-volume solvers for unstructured grids fitted to the particles is not straightforward, still has been tried by some researchers [9,10,11] and for dense packings, particularly of spheres, modifications of the original geometry (artificial gaps or bridges between particles) are necessary for generation of computational grids with sufficient quality [12,13]. Other popular approach is a simulation of flow with Lattice Boltzmann Method (LBM) [14,15]; see particularly [1,16] where LBM was applied to supply flow information for calculation of tortuosity of simple porous media.
In the present work, another method is proposed in which flow through a random packing is simulated with the classical second-order finite-volume solver for structured grids. The complex domain geometry is handled with the immersed boundary method originated by Peskin [17]. The paper is organised as follows. Section 2 contains the description of the employed methods -for numerical generation of the geometry of random packings (cylinders and rings), for flow calculations and estimation of the tortuosity parameter given the flow velocity field. In Section 3 sample results are presented and discussed. The paper closes with Summary in Section 4.

Numerical generation of a random packed bed geometry
The geometry of random packed beds can be obtained experimentally from real beds with the use of computer tomography techniques [18] but here a different approach is employed -Discrete Element Method (DEM) simulation using method presented in [19]. The main idea of the method is the simulation of the process in which the particles are sequentially dumped into a cylindrical container (column) and move due to gravity and reaction forces from other particles or the column wall. When The density of the packing can be increased by inclusion of additional force acting on the falling particles in the direction from the column's axis to its walls. Without this force the void fraction of the generated beds is relatively large, in fact, the largest possible to obtain in experiments. As demonstrated in [19,20] the algorithm is able to reproduce the structure of real beds, in particular void fraction distribution [18,21], with a sufficient accuracy.
For further calculations, at first, beds consisting of 2000 particles of diameter D p and height H p (in case of rings the inner diameter was D i ) were generated in a cylindrical container of diameter D c . As the base of the container introduces a distinct ordering of the particles [20], for the flow simulation only a part of the bed of height L was selected positioned over 5D p above the container's bottom to obtain more representative sample of the packing.

Flow solver and immersed boundary method
The fluid flow is solved with the in-house finite volume code developed at the Institute of Thermal Machinery based on regular MAC staggered grid and projection method [22,23,24]. Denoting flow velocity of the incompressible flow as [ , , ] x y z u u u = u , the flow equations take the form: where ρ -density of the fluid.
Time discretisation of equation (4) is performed with second order Runge-Kutta scheme. Spatial derivatives are also of second order -upwind scheme for the advection term and central scheme for the viscous term. The pressure equation in the projection method is solved with preconditioned conjugate gradient method (IPCG). More detailed description of the solver can be found in [25].
As calculations on regular structured grids are usually limited to simple domain geometries, the simulations in such complex domains as random packings require special treatment. In this case socalled immersed boundary method (IBM) is used [17,26,27,28]. The idea of the method is to take as the computational domain a regular area (e.g. a box) encompassing the actual flow domain and to modify equation (4) by inclusion additional source terms in order to enforce correct flow velocity on the "immersed" boundaries (here, on the particles' surface and the column wall). When the immersed boundaries are stationary with no-slip boundary condition, the flow velocity should be damped on them and it can be done after discretisation of the flow equation without explicit additional terms (direct forcing approach [26,28]). In the current solver, the damping of flow velocity is attained with the use of smooth damping function constructed from the actual geometry of the packing. As demonstrated in [25] this model leads to calculated values of pressure drop in packed beds of rings in a good agreement with empirical formulas. The flow configuration is shown in Figure 2

Passive advection of massless particles
When the steady state flow velocity u is known, the tortuosity of the streamlines can be found by passive advection of massless particles and calculation of the length of their trajectory (path). In this method the following equation is integrated: where x p is the position vector of a given particle. Initially the particles are distributed on a regular dense grid in the source plane located at distance L p over the top of the bed (see Figure 2). The equation (5)  the trajectories of different particles can significantly vary in length. Some of the paths end prematurely (see Figure 3a) because of the errors of integration or the fact the velocity field is extremely complex and not exactly divergence-free (Figure 3b). For the analysis of tortuosity only the trajectories covering at least 2/3 of the bed height were taken into account (about 80% of the total number of paths). The tortuosity calculated with this method will be further denoted as τ adv .

Method of mean velocities
Explicit calculation of trajectories is not the only way of estimation of bed tortuosity. In the literature another method was proposed in which tortuosity can be obtained directly from the flow velocity field. According to [16]: where the integration is performed over the void space in the packing (as the velocity equals 0 inside immersed boundaries integration over the whole computational domain gives the same result). This method will be further referred to as the method of mean velocities and the value estimated from formula (6) will be denoted as τ mvel .

Results and discussion
The methodology presented in Section 2 was used for calculation of tortuosities of four random packings: C0, C1, R0, R1 (see 2.1). The dimensions of the particles, the container and other characteristic lengths are shown in Table 1. In Table 2 the results were gathered for the tortuosities calculated with both methods: advection of particles (averaged value for sufficiently long paths) and mean velocities. We can see that for packings of lower density (C0 and R0) the agreement between the methods is very good. However, in the case of dense packings (C1 and R1) the discrepancies cannot be neglected. It must be mentioned that increased packing density is related not only to the reduced void fraction but also introduces additional ordering of the particles [18,20] which could influence the tortuosities of individual flow paths and the result of the averaging process, especially for the method based directly on the tortuosity definition (3).
In order to test the influence of grid density on the results, adv τ for the most demanding case R1 was calculated on four different grids: 80x112x80, 120x168x120, 160x192x160 and 200x280x160 giving 1.257, 1.291, 1.306 and 1.309 respectively. Good grid convergence justifies the use of the selected grid (160x192x160) also for other considered cases.  In the context of works connecting void fraction with tortuosity, it is worth noting that although void fraction of random packing of rings is prominently larger than that of packing of cylinders, the tortuosities are in fact comparable (see C0 and R0 cases). However, rings (and other hollow particles, in general) are definitely not typical particles considered in the search for ( ) τ τ ε = formula. As the obtained results show, the void space inside the rings contributes to the total void fraction of the bed but the complexity of flow paths is similar as in the case of full cylinders.
The method based on the advection of particles give information about tortuosity of each individual path, so beside the calculation of the averaged value representative to the whole bed more detailed characteristics, like the frequency distribution of particular tortuosity values, can also be determined ( Figure 4). In all analysed cases, the calculated values of τ adv are in major part in the range 1-1.6 (however, for C1 case values even as high as 6 can be found). Values in the range 1.2-1.25 are generally the most frequent. As could be expected, for the packings of lower density (C0, R0) the contribution from low values of τ adv (below the maximum) is more pronounced. This could result from wider pore channels extending over larger parts of the bed, in comparison with dense random packings.

Summary
The method proposed in this work can be employed for physically-base calculation of tortuosity of random packings using detailed information about flow velocity field in the low Reynolds number regime. The IBM method allows for application of a classical MAC solver to solution of the flow in complex geometries by simple modification of the projection method and generation of unstructured grid fitted to the boundaries is not necessary. Given steady state flow velocity field the tortuosity can be estimated directly from this field and for coarse packings of rings and cylinders the obtained values coincide with the ones determined explicitly from the definition of tortuosity. For dense packings, however, the results are markedly different.