## Abstract

We study the time and length scales of hydrodynamic dispersion in confined monodisperse sphere packings as a function of the conduit geometry. By a modified Jodrey–Tory algorithm, we generated packings at a bed porosity (interstitial void fraction) of *ε*=0.40 in conduits with circular, rectangular, or semicircular cross section of area 100*πd*^{2}_{p} (where *d*_{p} is the sphere diameter) and dimensions of about 20*d*_{p} (cylinder diameter) by 6553.6*d*_{p} (length), 25*d*_{p} by 12.5*d*_{p} (rectangle sides) by 8192*d*_{p} or 14.1*d*_{p} (radius of semicircle) by 8192*d*_{p}, respectively. The fluid-flow velocity field in the generated packings was calculated by the lattice Boltzmann method for Péclet numbers of up to 500, and convective–diffusive mass transport of 4×10^{6} inert tracers was modelled with a random-walk particle-tracking technique. We present lateral porosity and velocity distributions for all packings and monitor the time evolution of longitudinal dispersion up to the asymptotic (long-time) limit. The characteristic length scales for asymptotic behaviour are explained from the symmetry of each conduit’s velocity field. Finally, we quantify the influence of the confinement and of a specific conduit geometry on the velocity dependence of the asymptotic dispersion coefficients.

## 1. Introduction

Understanding transport phenomena in porous media, e.g. hydrodynamic dispersion, is important for many technological and environmental processes [1]. In particular, the topic of hydrodynamic dispersion in confined cylindrical packings of spherical particles at low cylinder-to-particle diameter ratios (≤20) has a long tradition in chemical engineering [2]. Under these conditions, longitudinal dispersion strongly depends on the flow maldistribution caused by the geometrical wall effect, which is characterized by damped oscillations of the interparticle porosity extending over a few particle diameters from the conduit’s inner surface towards the bulk of the packing [3]. This geometrical wall effect and the associated transcolumn velocity gradient result from the inability of spherical particles to pack closely against the rigid column wall. Thus, the presence of a confinement implies that at least two fundamental time and length scales exist in a sphere packing for longitudinal dispersion: the pore (short-time) scale and the confinement or transcolumn (long-time) scale [4,5].

The precise nature of how the microstructure of a packed bed and the geometry of the confinement affect flow heterogeneity, transverse equilibration and the resulting longitudinal dispersion is still largely unresolved. Specifically, packings in non-cylindrical conduits have rarely been studied [6,7], although their practical importance has steadily increased with the emerging lab-on-a-chip technology, and even though their unpacked counterparts, open channels of various cross-sectional geometries, have been extensively investigated by numerical analysis methods [8]. However, the results for open channels have only limited implication for confined packings, because here the packing microstructure determines the time and length scales that govern flow and dispersion.

Detailed three-dimensional numerical simulations of flow and transport in sphere packings are particularly suited to the challenge of resolving the central structure–transport relationships, because this approach allows us to systematically study relevant parameters, such as the shape and size of the particles, as well as conduit shape and packing length [4; 5,9–16]. Transient dispersion is recorded easily, quantifying time and length scales required for the attainment of asymptotic dispersion.

In this work, we report the results of a high-resolution numerical analysis of flow and mass transport in computer-generated packings of monosized, solid spheres confined in cylindrical as well as non-cylindrical conduits. The careful preparation of confined packings of sufficient length is indispensable to realize laterally fully equilibrated wall and corner effects. Our work quantifies the impact of the confinement by comparison with bulk packings (which mimic infinitely wide, random packings without confinement) and the influence of the confinement’s cross-sectional geometry on flow heterogeneity and the resulting transient and asymptotic dispersion as a function of the flow velocity.

## 2. Simulation methods

### (a) Packing generation

Packings of solid, impermeable, monosized spheres were generated by a modified Jodrey–Tory algorithm [12] at a bed porosity (interparticle porosity or average interstitial void fraction) of *ε*=0.40. Bulk (unconfined) packings of about 9.53*d*_{p}×9.53*d*_{p}×76.43*d*_{p} were periodic along all three directions. Confined packings (periodic in the longitudinal direction) were generated in containers of circular, rectangular (with a side aspect ratio of 2 : 1) or semicircular cross section of area 100*πd*^{2}_{p}. Packing dimensions were 20*d*_{p} (cylinder diameter) by 6553.6*d*_{p} (length), 25*d*_{p} by 12.5*d*_{p} (rectangle sides) by 8192*d*_{p} and 14.1*d*_{p} (radius of semicircle) by 8192*d*_{p}. Packings were discretized with a spatial resolution of 30 nodes per *d*_{p}, sufficient for the accurate simulation of fluid flow [10] and mass transport [13]. Discrete lattices had dimensions of 286×286×2293 nodes (bulk), 600×600×196 608 nodes (circular), 752×376×245 760 nodes (rectangular) and 848×424×245 760 nodes (semicircular). We used five random realizations of each packing type, resulting in a total of 20 generated packings. Results given for a packing type refer to the mean calculated from its five random realizations.

### (b) Simulation of fluid flow

Low Reynolds number flow (*Re*∼10^{−2}) of an incompressible fluid in the interparticle void space was simulated by the lattice Boltzmann method (LBM) [17]. We used the D3Q19 LBM model with a Bhatnagar–Gross–Krook (BGK) collision operator. Simulations were performed at very low lattice Mach numbers (*Ma*<10^{−3}) to minimize the compressibility error introduced by the LBM as *Ma* approaches *c*_{s}, the speed of sound in the system. A bounce-back (BB) rule was employed to implement the no-slip boundary condition at the solid–liquid interface (the surfaces of spheres and confinement). Because BB may introduce a significant error to pore-scale simulations if the value of the lattice viscosity *ν* differs from 1/6 [18], we used *ν*=1/6 in all our simulations. The selected approach (LBM–BGK with BB) and a grid resolution of 30 nodes per *d*_{p} allowed us to simulate hydrodynamic dispersion (using flow fields calculated by LBM) with a maximal error of 7 per cent in the asymptotic dispersion coefficients relative to a grid resolution of 90 nodes per *d*_{p}.

The D3Q19 lattice requires storage of 19 lattice links per lattice node, which resulted in computer memory requirements of *ca* 13 GB for bulk packings and of 5–6 TB for confined packings. The simulation of one velocity field required 1500 LBM iterations and took *ca* 0.8 h on 16 384 cores of a Blue Gene/P system (JUGENE at Forschungszentrum Jülich, Germany).

### (c) Simulation of mass transport

Advective–diffusive mass transport of conservative tracers in the interparticle void space was simulated using the random-walk particle-tracking technique [19]. The idea of the method is to distribute a large amount of tracers (*N*=4×10^{6} in our study) in the volume of interest and to track target quantities (tracer trajectories, statistical moments of the tracer ensemble, etc.) over time. In our simulations, tracer displacement during a time step *dt* is the result of two contributions, advective and diffusive. The former is derived from the flow velocity field (calculated by LBM); the latter is generated using normally distributed random numbers with zero mean and a variance of , where *D*_{m} is the tracer diffusion coefficient in the fluid. The hydrodynamic longitudinal dispersion coefficient was calculated as
2.1where Δ*r*_{zi} and 〈Δ*r*_{z}〉 are the corresponding Cartesian components of the longitudinal displacement of the *i*th tracer and the average longitudinal displacement of the tracer ensemble, respectively, after time *t*. Transverse dispersion coefficients *D*_{T}(*t*) were calculated similarly. The time step was chosen so as to restrict the maximal tracer displacement to *dl*/2, where *dl* is the lattice spacing. Mass transport simulations were performed on 8192 processor cores of the Blue Gene/P system and took *ca* 768 h for all confined packings.

Fluid flow and hydrodynamic dispersion were studied over a wide range of reduced (dimensionless) velocities or Péclet numbers, *Pe*=*u*_{av}*d*_{p}/*D*_{m} (where *u*_{av} is the average flow velocity), covering 0.1≤*Pe*≤500. The program realization of all algorithms was implemented in C language using the Message Passing Interface library. Further details on these in-house written codes, particularly their performance scaling, can be found elsewhere [20].

## 3. Results and discussion

### (a) Lateral porosity and velocity distributions

Figure 1*a* shows the front view of the confined sphere packings generated with a cross-sectional area of 100*πd*^{2}_{p} and a bed porosity of *ε*=0.40. The circular packing has a diameter of 20*d*_{p}, the side lengths of the rectangular packing are 25 and 12.5*d*_{p}, and the radius of the semicircular packing is 14.1*d*_{p}. The packing dimensions and porosity were chosen to reflect the conduit-to-particle-size ratios and relatively dense packings of the chromatographic beds used in capillary and microchip separations [21]. The lateral porosity distributions (figure 1*b*) reveal that particles are highly ordered in the near-wall region; the sphere centres of the first particle layer from the wall form a sharply defined line. Large, periodic fluctuations of local porosity are present over a distance of 3–4*d*_{p} from the wall, before the distribution approaches the randomness typical of the packing bulk. All confined packings display lateral porosity oscillations, but the amplitudes are larger for non-cylindrical conduits, because these contain corners and their symmetry is reduced compared with cylindrical conduits [6]. The velocity profiles in figure 1*c*, calculated for *Pe*=100, correlate well with the porosity profiles of figure 1*b*. Regions of higher- and lower-than-average velocity match the maxima and minima, respectively, of the porosity distributions. As anticipated from the impossibility to pack conduit corners tightly with monosized spheres, regions of increased velocity appear in the corners of the non-cylindrical packings.

### (b) Transient longitudinal dispersion

Transient behaviour of the longitudinal dispersion coefficients towards their asymptotic values was monitored for the confined packings, expecting a transcolumn contribution from the wall effects (cf. figure 1). Figure 2 shows the development of *D*_{L}(*t*)/*D*_{m} for *Pe*=100. Elapsed time in figure 2 has been normalized through the transverse dispersive time , where the value for *D*_{T}(*Pe*) was taken from the bulk packings and *L*_{C} is a characteristic transverse length for each conduit geometry. Based on the lateral porosity and velocity distributions of a particular conduit geometry, this length is the straight distance through a packing cross section that needs to be covered by a tracer for equilibration, i.e. to experience all the different velocities in the velocity field. Under explicit consideration of a velocity field’s symmetry, the characteristic length is the maximum distance between different velocities.

The highest velocities are found along the wall of cylindrical conduits and in the corners of non-cylindrical conduits. Complete exchange between different velocities in the fluid is achieved by covering the lateral distance from wall to centre in a circular packing, from corner to centre in a rectangular packing and from corner to apex in a semicircular packing. Thus, the characteristic transverse length on the macroscopic (conduit cross-sectional) scale is the radius for circular conduits, the half-diagonal for rectangular conduits and the distance from the apex to a corner for semicircular conduits.

Our definitions of the characteristic transverse lengths *L*_{C} for the confined packings are validated by the longitudinal dispersion data (figure 2), which demonstrate attainment of asymptotic behaviour at *t*_{D}≈1 for all conduit geometries. This dispersive time of one is the time span after which tracers are laterally dispersed by one *L*_{C}. For the cylindrical packing, for example, this result confirms that the macroscopic flow heterogeneity in a monodisperse sphere packing caused by the cylindrical confinement adds a transcolumn contribution to the dispersion, which—owing to cylindrical symmetry—requires a lateral equilibration on the scale of the cylinder radius (*L*_{C}=*d*_{c}/2, where *d*_{c} is the cylinder diameter). The observed asymptotic time scale of *t*≈(*d*_{c}/2)^{2}/2*D*_{T}(*Pe*) for longitudinal dispersion is reminiscent of the classical Taylor–Aris dispersion in laminar (Poiseuille) flow through an empty cylindrical column (open tube), where the asymptotic time scale is proportional to *t*≈(*d*_{c}/2)^{2}/2*D*_{m} [1]. The decreasing symmetry from circular to rectangular to semicircular cross section of the conduits (from bottom to top in figure 2) results in an increased characteristic transverse length of the velocity field at equal cross-sectional area.

It should be noted that experimental studies are often carried out under non-equilibrium conditions, where solute residence times in a packing are insufficient to fully relax the transcolumn contributions to hydrodynamic dispersion. Thus, if the complete process shown in figure 2 is truncated (i.e. asymptotic longitudinal dispersion is not attained), it may be incorrectly concluded that wall effects are absent. Figure 2 demonstrates that, owing to the transcolumn velocity bias, which depends on the specific conduit geometry, the observation of truly asymptotic dispersion behaviour and its velocity dependence in confined systems require sufficiently long simulation times (and packings).

### (c) Velocity dependence of asymptotic dispersion

Figure 3 compares the velocity dependence of the asymptotic longitudinal dispersion coefficients (for , cf. figure 2) of the confined packings with the corresponding behaviour of bulk (unconfined) packings. In hydrodynamics, the velocity dependence of the dispersion coefficient is usually described by a power law [1]. For *Pe*≥10, we neglect in a first approximation the effect of molecular diffusion on longitudinal dispersion and write for the power law in the convection-dominated regime (10≤*Pe*≤500): [9]. From a least-squares fit (*R*^{2}>0.9999) to the datasets in figure 3 we obtained the values for the parameters *A* and *α* for each packing type.

The exponent for the bulk packing (*α*=1.191) agrees excellently with experimental data (*α*∼1.2; [1]) for various unconsolidated media (bead packs, sand packs, sandstones). For confined packings, the exponents are decidedly larger, although they remain similar to each other. The exponent for the confined cylindrical packing (*α*=1.333) is close to the experimental value of *α*=1.37 reported by Seymour & Callaghan [22] for a similar cylinder-to-particle diameter ratio (approx. 20), though at a slightly higher bed porosity (*ε*=0.44). The coefficient *A* depends on the heterogeneity of the pore space [1]. We observe (i) a large increase in *A* from bulk to confined packings, and (ii) a systematic increase in *A* for the confined packings with decreasing conduit symmetry (i.e. increasing lateral equilibration distance for velocity heterogeneities). This behaviour agrees with the experimental trend observed by Jung *et al.* [21] for longitudinal dispersion in particulate microchip packings contained in channels of different cross-sectional geometries, where dispersion was found to increase with decreasing conduit symmetry (at similar bed porosity).

The quantitative analysis in figure 3 confirms the impact of the confinement-inflicted transcolumn velocity bias on dispersion as well as the role of the particular conduit shape (decreasing symmetry from circular to rectangular to semicircular). Most important, the dispersion coefficients of the confined packings are up to 10 times larger than those of the bulk packings.

## 4. Conclusions

The presented numerical approach allows the large-scale simulation of flow and dispersion in three-dimensional confined packings in containers of arbitrary shape. It provides quantitative information on the spatial distribution of simulated velocities and on the resulting transient as well as asymptotic dispersion. We used this feature to analyse the impact of confinement and of the conduit shape on the velocity profiles and hydrodynamic dispersion in a pressure-driven flow through monodisperse sphere packings over a wide range of velocities. Our results on the velocity dependence of the asymptotic longitudinal dispersion coefficient show that the change from bulk to confined packings increases dispersion substantially owing to the transcolumn velocity bias caused by the wall effect. Our data also reveal two important aspects about the conduit geometry: (i) channels of advanced flow velocity reside in the corners of the non-cylindrical conduits, and (ii) the reduced symmetry of non-cylindrical conduits translates to a longer characteristic length of the tracer for lateral equilibration between different velocities; consequently, hydrodynamic dispersion in non-cylindrical packings is larger than in cylindrical packings of equal cross-sectional area. Among the investigated conduit shapes, the semicircular packing has the worst dispersion characteristics, because it has the lowest symmetry and a decidedly longer characteristic transverse length.

## Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft DFG (Bonn, Germany) under grant TA 268/4-1. We thank the DEISA Consortium (www.deisa.eu), co-funded through the EU FP6 project RI-031513 and the FP7 project RI-222919, for support within the DEISA Extreme Computing Initiative.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society