The macroscopic properties of two-dimensional random periodic packs of polydisperse cylinders are investigated by means of numerical methods. We solve the unsteady, two-dimensional Navier-Stokes equations on a staggered Cartesian grid and use the immersed boundary method to treat internal flow boundaries. The effects of porosity, polydispersity, and Reynolds numbers on the macroscopic permeability are studied. For small Reynolds numbers, we show that the permeability can be correlated to the underlying microstructure by means of a suitably defined statistical descriptor, the mean shortest Delaunay edge. With proper scaling, the results for polydisperse cylinders collapse onto the data for monodisperse cylinders, which can then be fitted with a universal curve. We also carry out a statistical analysis of the permeability computed for 500 samples and show that rare events, where the permeability lies outside the mean plus/minus three times the standard deviation, are possible. Finally, for larger Reynolds numbers, we show that a modified Forchheimer equation can characterize the flow.

## I. INTRODUCTION

Fluid flow in porous media has been studied for centuries, yet due to wide variations in typical morphologies, it is still difficult to predict flow properties in many materials. In the past, correlations have long been popular for predicting the pressure drop in flows through porous media. For example, Darcy^{1} was one of the earlier pioneers in the subject, who carried out extensive experiments and proposed a linear correlation relating the superficial fluid velocity to the pressure drop across the media using the proportionality *μ*/*K*. Here, *μ* is the fluid viscosity and *K* is the permeability. The permeability *K* is found by experiment or from another correlation and is a material property of porous media. Darcy’s law is valid at the macroscale and defines a relationship between the macroscopic velocity and the pressure; the term macroscale refers to spatial averaging on a scale large compared to the pore geometry. Darcy’s law is valid as long as the inertial forces are negligible compared to viscous forces. In the inertial regime, other correlations are popular such as Ergun’s equation,^{2} a generalization of the Forchheimer equation.^{3}

Although Darcy’s law in the non-inertial regime appears straightforward to use in predicting macroscopic flow quantities, its use is complicated by the fact that typical porous media has complex microstructures. Little is known about how the microstructure (polydispersity; particle shape; porosity; inter-particle spacing; higher-order statistics) affects macroscopic quantities such as fluid flow and permeability.

Numerical simulations are becoming popular as an alternative to experiments. For a precise description of the fluid flow through a porous medium, the Navier-Stokes equations are solved in a predefined pore geometry. However, because of computational limitations, most studies are limited to ordered packs of mono- or bidisperse cylinders/disks (2-D) and spheres (3-D) in square, triangular, or hexagonal arrays. For random packs of particles, little is known about how the microstructure affects permeability. There have been limited studies on the effect of polydispersity and particle shape on permeability, and we review some of these below.

One example of complex microstructures is natural sandstone. The work presented in Ref. 4 investigate three-dimensional pore space geometry for three natural sandstones. The geometries were reconstructed using microtomography. Using local porosity distributions and local percolation probabilities, the authors performed quantitative comparisons between experiments with microscopic pore space models across the three samples. However, due to the complexity of the microstructure, fluid flow through the geometries was not considered.

Another study that used microtomography to generate the microstructure is that of Ref. 5. The authors carried out a comparison of digitized pore structures (microtomography images of real lithofacies samples) and numerical reconstruction of rocks to predict material transport properties of complex, heterogeneous rock. They solve Stokes flow in the pore region using a lattice Boltzmann algorithm and compute, among other quantities, the permeability. The results show good agreement between the tomography images and reconstructed samples, demonstrating the potential use of digitized pore samples and numerical simulations.

Still another study is that of Ref. 6, who used both microtomography of a sand sample and an artificial method using a Boolean model for generating various microstructures. The Boolean model is based on a random distribution of overlapping ellipsoids whose size and shape were chosen according to some criteria. The work mainly focused on geometric parameters called Minkowski functionals (volume, surface, curvature, and connectivity) and their effects on permeability. The flow through the samples was solved using a lattice-Boltzmann method. The results show that the Minkowski functionals are not sufficient for the characterization of the geometric properties, whereas the permeability of the scanned and artificial samples were in good agreement.

Another example is that of Ref. 7, who used numerical simulations to show that the pressure drop across a bed of coal particles (randomly shaped polyhedra generated using Voronoi tessellation) changes by more than a factor of two when the equivalent particle diameter is kept constant while the polydispersivity of the pack changes. However, this work was limited to only two particle size distributions.

Garcia *et al.*^{8} computed the permeability as a function of particle shape and polydispersity for five three-dimensional random packs of non-spherical particles. The particle shape and grain size distribution for three samples were chosen to mimic real sand, while monodisperse and bidisperse packs were generated for comparison. Stokes flow was solved in the three-dimensional pore region by use of a finite element method. Their results show that the permeability is strongly affected by porosity and mean grain size, but particle shape and size distributions have only a small effect.

The permeability and tortuosity were computed for random packs of polydisperse disks in Ref. 9. The particle size distributions were calculated by a Weibull distribution, and random packs were then generated by a discrete element method. The number of different particle sizes used in the study was chosen according to a gradation curve (cumulative distribution function); only four particle diameters were considered for the steepest gradation, while six particle diameters were taken for the slowest gradation curve. The fluid flow in the pore regions was solved by the lattice Boltzmann method, and the results were correlated via either the Hazen or the Carman-Kozeny equation. The results demonstrated the effects of porosity and specific surface ratio on the permeability.

Despite the pioneering works mentioned above, there are a number of shortcomings. The most obvious one is the lack of defining a proper representative volume element (RVE). As pointed out in Ref. 10, when the RVE is sufficiently large or under the assumption of ergodicity, when a sufficiently large ensemble size is used for relatively smaller RVEs, the mean permeability should converge. Thus, using a single realization of a random complex microstructure does not allow one to determine the mean and standard deviations of global quantities such as permeability with any degree of confidence. In addition, using realizations (either microtomography or computer generated) on domains with solid walls rather than periodic boundaries can lead to boundary effects, which can have an order one effect on permeability. To address these issues, we use ensemble averaging on random packs of polydisperse cylinders (disks in two dimensions, or flow across infinitely long fibers in fibrous materials) with periodic boundaries to compute the permeability.

An interesting concept in coupling the microscale to the macroscale is to correlate the nondimensional permeability to some statistical descriptor of the underlying microstructure. This was recently considered for two-dimensional random packs of monodisperse cylinders.^{10-15} In the first work, the authors find that a correlation exists between the permeability and the mean nearest inter-fiber spacing, which takes into account the fluid space between two fibers. In the second work, the authors characterize the flow using Delaunay triangulation. In the latter works, the authors characterize the packs using particle nearest neighbor distances, Delaunay triangulation, and Delaunay edges. Physically, the line segments formed by nearest neighbors or shortest Delaunay edges define, roughly speaking, the edges of two-dimensional channels where the fluid will most likely flow (fast). By examining random packs, the authors show that (i) the permeability does not scale with the first nearest neighbor length scale; (ii) the second nearest neighbor is roughly equal to the shortest Delaunay edges; and (iii) a universal curve exists that relates the permeability to the second nearest neighbor.

In this paper, we solve for fluid flow in random periodic packs of polydisperse cylinders. We solve the unsteady Navier-Stokes equations on a Cartesian grid and use the immersed boundary method to treat internal flow boundaries. The immersed boundary method imposes the boundary conditions on the body surface by implicitly updating the flow variables on ghost cells which are located just inside the body. Since this method does not require body-conforming grids, it allows the use of regular Cartesian grids, removing the time-consuming process of building body-conforming grids. It is advantageous especially when there is a need to handle complex, multiple shapes in the flow and/or moving boundaries; although moving boundaries are not considered in the current work. For a recent application of the immersed boundary method to in-line and staggered arrangement of infinitely extending square rods, see Refs. 16–18. Although three-dimensional simulations can be carried out on today’s computers, we consider only two-dimensional packs. The reason for this is (i) to see if the immersed boundary method is a viable alternative to body-conforming grids; (ii) to get a quicker turn around in terms of wall clock time to see if using a statistical descriptor, such as Delaunay edges, would be a fruitful path forward in three-dimensions; and (iii) two-dimensional studies are of interest in and of themselves, for example, for certain fibrous materials where the flow is across the fibers. In future work, we plan on carrying out three-dimensional simulations for random packs of polydisperse particles.

The structure of the paper is as follows. In Sec. II, we present the governing equations and the numerical method. Results are presented in Sec. III where we discuss our packing algorithm used to generate random packs, the effect of porosity and polydispersity on permeability, and present a statistical analysis of 500 packs. We also present results at higher Reynolds numbers. Finally, conclusions are given in Sec. IV.

## II. PROBLEM FORMULATION AND NUMERICAL METHOD

We solve by means of numerical methods the unsteady incompressible Navier-Stokes equations in two dimensions

where *i*, *j* = 1, 2, (*u*_{1}, *u*_{2}) = (*u*, *v*) are the velocity components in the *x*, *y*-directions, respectively, *p* the pressure, and *Re* the Reynolds number. The variables are scaled by a reference velocity *U*_{0}, lengths by *L _{p}*, and pressure by $\rho U 0 2 $. The Reynolds number is then defined as

*Re*=

*ρU*

_{0}

*L*/

_{p}*μ*, where

*μ*is the fluid viscosity assumed constant,

*ρ*the fluid density, and

*L*a length scale associated with the microstructure, such as a particle diameter.

_{p}For the temporal integration of (1), the fractional-step method with the low-dissipation and low-dispersion third-order Runge-Kutta scheme of Ref. 19 is employed. The convective terms in (1) are discretized by a 3rd-order upwinding WENO scheme,^{20} the diffusion terms by a 4th-order central difference scheme, and the pressure term by a 2nd-order central difference. At each stage, the Poisson equation for the pressure is solved by a standard iterative multigrid solver.

The computational domain is discretized using a staggered Cartesian grid with uniform mesh size. The staggered grid is imposed to prevent odd-even decoupling of the pressure. The *u*-velocity is located at the face-center on horizontal sides, the *v*-velocity is on vertical sides, and the pressure is placed at corners of each cell; Figure 1 shows the location of the velocity components and the pressure. The latter is done to take advantage of standard multigrid solvers where the boundary conditions are applied along the unit cell edges and thus no extrapolation of the pressure at the boundaries is required.

Internal flow boundaries in the flow field are treated using the immersed boundary method.^{21,22} The term “immersed boundary method” refers to any method that solves the Navier-Stokes equations with embedded boundaries on grids that do not conform to the shape of these boundaries. The major advantage of using immersed boundary methods over conventional body-fitted (conformal) methods is the ease with which arbitrary-shaped immersed bodies can be treated; using body-fitted meshes for complicated geometries can quickly lead to deterioration in grid quality and can have a negative impact on accuracy and convergence of the results.

Implementation of the immersed boundary method begins with identifying the fluid, solid, and ghost nodes (GNs). The identification of a fluid or solid node begins by defining a level set function *f*(*x*, *y*) such that a grid point lies in the fluid domain if *f* > 0 and in the solid domain if *f* < 0. The fluid/solid boundaries are therefore defined by the level curve *f* = 0. To adequately represent the flow boundaries and to interpolate ghost nodes across the fluid/solid boundaries (to be discussed below), the level set function needs to reside on a grid that is much finer than the flow domain grid. Numerical tests suggest that a grid that is four times finer than the flow domain is sufficient for convergence; in all of our simulations, we take a grid that is six times finer than the flow domain grid. A fluid/solid node is a node that lies in the fluid/solid region, respectively. A GN is defined as a node that is located inside the solid body, and at least one of the adjacent north, south, east, and west nodes is in the fluid. Since each flow variable is located on a different node due to the staggered grids, the property of the node (fluid, solid, or ghost) is determined only by its coordinates. An image point (IP) is defined on the line segment extended from each ghost node normal to the body surface and lies in the fluid domain. A cross point of the line and the body surface is defined as a boundary intercept (BI) located in the midway between IP and GN. A flow variable on IP is calculated by a bilinear interpolation from the values on the surrounding four nodes. A linear interpolation is then applied to determine the values on GN from those on IP and BI (i.e., the given boundary conditions). At each sub-step of the temporal integration, the variables at the fluid nodes are updated by solving the governing equations and those at the ghost nodes are then updated by the above procedure when the boundary conditions are applied. The variables at the solid nodes are simply set to zero. Figure 1 shows a schematic of the ghost nodes, the image points, and the boundary intercepts. See Refs. 21 and 22 for more details.

Extensive validation and verification (V&V) tests of the numerical method were carried out and reported in Refs. 23 and 24; here, we summarize the findings for completeness. To test the solver without immersed boundaries, we solve the Taylor-Green vortex problem. Numerical solutions are compared to the exact analytical solution, and results show convergence is second-order for velocity and pressure, as expected. In addition to the Taylor-Green problem, we also solve the driven cavity problem. Numerical results are compared with the numerical results from other publications. The results show excellent comparison with differences of less than 1% for Reynolds numbers up to 3200. The coupled Navier-Stokes solver and immersed boundary method was then tested by carrying out simulations of flow past a cylinder and a square rod, which are standard test problems found in the literature. For Reynolds numbers of 20 and 40, the flow is steady, and the computed drag coefficient and wake bubble length compare to within 2% of that computed numerically by others. For Reynolds number of 80, the flow is unsteady, and the computed drag and Strouhal number compare well with other published results. To test the solver in computing the permeability, a set of V&V studies was carried out for flow through monodisperse and bidisperse packs of cylinders in square and hexagonal arrays. Using a modified Forchheimer equation, the friction factor, when plotted against a modified Reynolds number, was seen to collapse onto a single curve for Reynolds numbers up to about 2. These results were found to be in excellent agreement with those of Ref. 25. Additional comparisons are made for random packs of monodispersed cylinders. The computed permeability as a function of porosity was compared to the results of Refs. 10 and 26, and again excellent agreement was obtained. The interested reader should consult Refs. 23 and 24 for further details.

## III. FLOW THROUGH RANDOM PACKS OF POLYDISPERSE CYLINDERS

In this section, we present results for fluid flow through two-dimensional random packs of polydisperse cylinders. To generate the packs, we use our in-house packing code, which is briefly described below. For a given microstructure, we solve the two-dimensional Navier-Stokes equations (1). The velocities are assumed to be periodic in both *x*- and *y*-directions, while the pressure in the *y*-direction is periodic. The flow is driven by applying a pressure drop between the inlet and outlet as a Dirichlet boundary condition. The domain size is determined from the porosity. The flow is solved numerically until it reaches steady state. The grid resolution used in this section is *dx* = *dy* ∼ [0.0065 − 0.01], depending on the value of the porosity. These resolutions are selected so that roughly a two times finer resolution yields differences of less that 1.0% (results not shown).

### A. Packing algorithm

In previous work, we described a packing algorithm that can generate packs of polydisperse spheres for mesoscale simulations of solid propellants in an efficient manner.^{27-29} The packing code begins with the basic Lubachevsky-Stillinger (LS) dynamic packing algorithm.^{30} The LS algorithm starts by placing *N* spheres with zero radii (points) at random locations inside a domain of interest. When the spheres have dimensional units associated with the diameters, the packing algorithm uses the largest diameter for nondimensionalization. The packing domain is then scaled to be on the unit cube. The spheres are given random velocities sampled from a Maxwell-Boltzmann distribution at a fixed temperature and are allowed to grow at a specified growth rate. The particles undergo classical (super-)elastic collision dynamics as they grow to fill the space in the domain. The algorithm stops when either a specified packing fraction is reached or when a specified jamming criterion is met. The algorithm was improved by incorporating an event-driven molecular dynamics approach and a hierarchical cell structure to reduce the overall computational time.^{29} Validation of the packing code against existing experimental data of spheres was also reported in Ref. 29. The packing algorithm was recently extended to pack non-spherical convex shapes using a novel method based on level sets.^{31,32} This algorithm extension was necessary because most particles of heterogeneous materials often have non-spherical shapes.

Here, we use our packing code to generate two-dimensional random periodic packs of polydisperse cylinders at various porosities. We specify the minimum distance between cylinders to be 0.1 of the radius of the largest particle, necessary for two-dimensional packs to avoid complete blockage of the flow.^{10}

For each pack considered, we compute the permeability. The permeability *K* is scaled by the square of an appropriate length scale *L _{p}* and is determined using Darcy’s law. That is, for a given pressure drop across the unit cell, the superficial (or average) velocity is first determined from the numerical simulation, and then the nondimensional permeability

*K*′ is solved for using Darcy’s law

where *Re* = *ρU*_{0}*L _{p}*/

*μ*is the particle Reynolds number with fluid density

*ρ*, velocity scale

*U*

_{0}= 〈

*u*〉, and viscosity

*μ*; 〈

*u*〉 is the non-dimensional superficial velocity; Δ

_{s}*p*is the non-dimensional pressure drop imposed across the unit cell in the transverse

*x*-direction; and

*L*=

*L*

^{∗}/

*L*is the nondimensional length, with

_{p}*L*

^{∗}the dimensional length of the pack domain in the flow direction.

*L*is taken to be the particle diameter

_{p}*D*for monodisperse cylinders, while for polydisperse cylinders it is taken to be

Here, *M* is the total number of unique particle sizes and *n _{i}* is the number fraction for particle

*i*with diameter

*D*. The total number of cylinders in a pack is denoted by

_{i}*N*.

_{p}### B. Representative volume element

When two- and three-dimensional random packs are used in numerical simulations, it is necessary to carry out ensemble averaging to account for the effect of orientations, spacing, etc., on quantities of interest such as permeability. The ensemble average of the permeability and its standard deviation are defined as

where *N _{r}* is the ensemble size and

*K*is the computed permeability for the

_{i}*i*th pack. The corresponding mean nondimensional permeability is therefore defined to be $ K \u2032 =\u3008K\u3009/ L p 2 $. As pointed out in Ref. 10, when RVE is sufficiently large or under the assumption of ergodicity, when a sufficiently large ensemble size is used for relatively smaller RVEs, the mean permeability should converge. Note that the larger the value of

*N*, the larger the domain for a given porosity. In all of the results presented below, we fix the total number of cylinders to be

_{p}*N*= 196, and we fix the ensemble size to be

_{p}*N*= 10. These values are based on previous work for random packs of monodisperse cylinders

_{r}^{24}and are deemed sufficient for our purposes.

### C. Pack properties

The particle size distribution is generated using a Gaussian distribution, defined by

where *f*(*D*) is the probability with standard deviation *D*_{σ} and mean $ D \u0304 $. The mean diameter is scaled such that $ D \u0304 =1.0$ and we use nine diameters; these are chosen to be 1.0, 1.0 ± 0.5*D*_{σ}, 1.0 ± 1.0*D*_{σ}, 1.0 ± 1.5*D*_{σ}, and 1.0 ± 2.0*D*_{σ}. To investigate the effect of polydispersity, we take *D*_{σ} = 0, 0.125, 0.25, and 0.375. The case of *D*_{σ} = 0.375 corresponds to a large-to-small diameter ratio of 7, a value typical of most applications. Thus, in this study, we do not investigate highly polydisperse distributions. Note that *D*_{σ} = 0.0 indicates a monodisperse pack. For each value of *D*_{σ}, packs with porosity of *σ* = 0.45, 0.60, and 0.70 are generated. The distributions used in this study are plotted in Figure 2(a). The corresponding diameters and number of cylinders for each value of *D*_{σ} are tabulated in Table I. A representative random pack is shown in Figure 2(b) for *D*_{σ} = 0.25 and porosity *σ* = 0.45.

Diameter . | . | ||
---|---|---|---|

D_{σ} = 0.125
. | D_{σ} = 0.25
. | D_{σ} = 0.375
. | No. of cylinders . |

0.25 | 0.50 | 0.75 | 5 |

0.44 | 0.63 | 0.81 | 13 |

0.63 | 0.75 | 0.88 | 24 |

0.81 | 0.88 | 0.94 | 36 |

1.00 | 1.00 | 1.00 | 40 |

1.19 | 1.13 | 1.06 | 36 |

1.38 | 1.25 | 1.13 | 24 |

1.56 | 1.38 | 1.19 | 13 |

1.75 | 1.50 | 1.25 | 5 |

Diameter . | . | ||
---|---|---|---|

D_{σ} = 0.125
. | D_{σ} = 0.25
. | D_{σ} = 0.375
. | No. of cylinders . |

0.25 | 0.50 | 0.75 | 5 |

0.44 | 0.63 | 0.81 | 13 |

0.63 | 0.75 | 0.88 | 24 |

0.81 | 0.88 | 0.94 | 36 |

1.00 | 1.00 | 1.00 | 40 |

1.19 | 1.13 | 1.06 | 36 |

1.38 | 1.25 | 1.13 | 24 |

1.56 | 1.38 | 1.19 | 13 |

1.75 | 1.50 | 1.25 | 5 |

### D. Permeability

The mean nondimensional permeability as a function of porosity and for various values of *D*_{σ} is plotted in Figure 3(a). The figure shows that the permeability increases as the porosity and polydispersity (variations in diameter) increase. From the figure, it is tempting to see if a scaling exists that collapses the data onto a single curve. And so we show in Figure 3(b), the results with the following scaling:

where $ K s \u2032 $ is the scaled nondimensional permeability and *s*_{1} = *s*_{1}(*D*_{σ}) is a quadratic function in *D*_{σ}. A least squares fit results in the values (*a*_{0}, *a*_{1}, *a*_{2}) = (1.016, − 0.489, − 0.595). With this scaling, the data essentially collapse onto the monodisperse results *D*_{σ} = 0.0.

In previous work, it was shown that the permeability can be related to the microstructure by means of some suitably define statistical descriptor. This was recently considered for two-dimensional random packs of monodisperse^{13,15} and bidisperse^{24} cylinders. In these works, the authors characterized the packs using Delaunay edges, defined as

Here, $ e 1 T $ represents the shortest edge of each Delaunay triangle, $ R \u02c6 $ represents the sum of the two radii that make up the shortest edge (the difference between $ e 1 T $ and $ R \u02c6 $ therefore represents the interparticle spacing), and *L _{p}* is given by Eq. (3). Physically, the line segments formed by the shortest Delaunay edges define, roughly speaking, the edges of two-dimensional channels where the fluid will most likely flow (fast). This can be better visualized by examining Figure 4 where contours of the velocity magnitude are plotted along with the mean shortest Delaunay edges.

The parameter $ \gamma 1 T $ as a function of porosity and for various values of *D*_{σ} is plotted in Figure 5(a). From the figure, we see that the value of $ \gamma 1 T $ increases as the porosity and polydispersity increase. Again, it is tempting to determine a scaling that will collapse the data onto the monodisperse case. We find that indeed such a scaling exists and is given by

with (*b*_{0}, *b*_{1}, *b*_{2}) = (1.006, − 0.144, − 0.824). Figure 5(b) plots the scaled results.

Finally, the mean nondimensional permeability as a function of $ \gamma 1 T $ and for various values of *D*_{σ} and *σ* is plotted in Figure 6(a) and the scaled values in Figure 6(b). In each figure, the prediction given by

is plotted as the solid line. This equation was proposed in Ref. 24 as a universal fit for random packs of monodisperse and bidisperse cylinders. The exponent of 2.5 in (9) is from lubrication theory for square and hexagonal arrays of monodisperse and bidisperse cylinders,^{33,34} and the exponent *ν* takes into account disordered arrays. Currently, there is no theory to determine *ν* for random packs of particles. It can be seen from Figure 6(b) that the scaled permeability is well predicted by (9).

### E. Statistical analysis and rare events

It is well known that when the number of random samples is sufficiently large, there may appear a case when the computed permeability is drastically off of the mean. Such a case is called a rare event, usually defined when a value is larger or smaller than the mean plus/minus three times the standard deviation.

To investigate if such events exists for random packs of polydisperse cylinders, the permeabilities of 500 samples with *σ* = 0.45 and *D*_{σ} = 0.25 are computed. The results are plotted along with the mean and standard deviations in Figure 7(a) for the nondimensional permeability and Figure 7(b) for $ \gamma 1 T $. The means and standard deviations are computed from Eq. (4). From the figure, we see that the permeability and $ \gamma 1 T $ at six points are found to be rare events. This suggests that rare events can exist for random packs and that rare events in the permeability do not necessarily correspond to that of $ \gamma 1 T $. We comment that a grid resolution was carried out for the rare event cases to ensure that the value of the permeability did not change with resolution; we find that doubling the grid in both directions only changed the value of the permeability by less than 2%. We also comment that we do not believe that the separation distance, used in our packing algorithm to avoid complete flow blockage, has any effect on the presence of rare events. We will revisit this issue in future work on three-dimensional packs, where a separation distance is not needed.

To analyze the microstructural behavior of rare events, pressure contours with permeability close to the mean (*K*′ = 0.00141) and close to a rare event (*K*′ = 0.00081) are shown in Figure 8. The mean Delaunay shortest edges are also plotted as white line segments. Note that the high pressure region (yellow) is narrower in the flow direction when the nondimensional permeability is close to the mean when compared to that of the rare event. This is probably due to the fact that the cylinders are more homogeneously packed when the permeability is close to the mean, whereas for the rare event, the pack has a large voidage region at roughly (*x*, *y*) = (−3, 6) and a dense region at roughly (*x*, *y*) = (7.5, − 8). The corresponding contour plot of the velocity magnitude (not shown) for the rare event shows a peak magnitude of roughly 0.07 in the dense region; this value is about twice that of the case where the permeability is close to the mean (cf. Figure 4). The corresponding streamlines are plotted in Figure 9.

The pressure distribution is also analyzed by plotting the mean pressure profile along the flow direction. The mean pressure in the fluid region is averaged over the *y*-axis as

where *L _{fluid}*(

*x*) is the vertical length of the fluid region at each

*x*-location, and

*ymin*and

*ymax*represent the height of the domain. A similar definition is used for the mean velocity magnitude $ u \xaf ( x ) $. Figure 10 plots the mean pressure $ p \xaf ( x ) $ and mean velocity magnitude $ u \xaf ( x ) $ as a function of

*x*with permeability close to the mean (

*K*′ = 0.00141), two rare events with permeability below −3

*σ*(

_{std}*K*′ = 0.00081 and

*K*′ = 0.00095), and a rare event with permeability greater than +3

*σ*(

_{std}*K*′ = 0.00183). Figure 10(a) shows that the mean pressure profile decreases roughly uniformly in

*x*. The profile with the high permeability does not show a distinct difference from the mean permeability profile. On the other hand, the mean pressure profiles corresponding to the two rare events with low permeability show regions of where the pressure drops rapidly as well as where the pressure remains relatively flat. This behavior in the pressure agrees with the above observation concerning the pressure contour. Figure 10(b) shows that the velocity magnitudes for the two rare events with low permeability lie about 40% below that for the mean permeability case.

Finally, to examine the statistical nature of the permeability, we show in Figure 11(a), the probability using a normal distribution, and in Figure 11(b), the probability using a Weibull distribution. For the normal distribution, the mean and standard deviations are found to be *μ* = 0.00142 and *σ _{std}* = 0.00012, respectively; for the Weibull distribution, the scale and shape parameters are found to be

*a*= 0.00147 and

*b*= 12.21, respectively. From the figures, we see that although the data are well represented by a normal distribution, the fit to the Weibull distribution is better, especially in the tails. Figure 12 compares the normal distribution (solid) to the Weibull distribution (dashed). Note that since the value of

*b*> 3.7, the Weibull distribution is negatively skewed (left tail). Therefore, it might be better to describe the data using a Weibull distribution rather than assuming a normal distribution.

### F. Effect of Reynolds number

It is well known that Darcy’s law (2) is only valid in the Stokes limit. At higher Reynolds numbers, the relation between the pressure drop ∇*p* and the superficial velocity 〈*u*〉 is found to be nonlinear due to inertial effects. To correct for this, Forchheimer proposed a quadratic correlation in the superficial velocity,^{3} given by

where *ρ* is the density of the fluid, *μ* the viscosity, *L* a length scale associated with the microstructure, *K* the permeability, and *β* an empirical, material-dependent parameter. The Forchheimer equation is written as the sum of a viscous term giving rise to the linear term and an inertial resistance giving the quadratic term. The Forchheimer equation is valid over a wider range of Reynolds number when compared to Darcy’s law (*β* = 0) but it requires experiments to determine the parameters. Ergun^{2} established a semi-empirical, generalized form of the Forchheimer equation for homogeneous bed of spheres. It is common to write Ergun’s (or equivalently Forchheimer’s) equation in terms of a friction factor and a Reynolds number. For example, Ergun’s equation becomes

where *D* is the mean particle diameter. The most relevant form of Forchheimer’s equation for our discussion is due to Papathanasiou *et al.*,^{25} given by

where $ f K $ is the modified friction factor, $R e K $ is the modified Reynolds number, and *F* is an empirical parameter. Based on square and hexagonal packings of disks (two-dimensional fibers), the authors showed that the best fit yields *F* = 0.08. Other similar definitions can be found in Table II of Ref. 14.

To determine the effect of Reynolds number on the permeability in the steady (viscous and inertial) regime, we compute the friction factor for random packs of polydisperse cylinders. We first show results using Ergun’s equation. Figure 13 plots the friction factor *f _{D}* as a function of Reynolds number

*Re*for various values of the porosity

_{p}*σ*and polydispersity parameter

*D*

_{σ}. Also shown in the figure is Ergun’s equations (12), with ±30% variations. We see that at the lower porosities of 0.45 and 0.6, Ergun’s equation under-predicts the numerical values by approximately 30% over the range of Reynolds numbers considered. At the higher porosity

*σ*= 0.7, the agreement between the numerical solutions and Ergun’s equation is significantly improved. This observation is consistent with that of monodisperse cylinders,

^{24}suggesting that Ergun’s equation is not a universal curve. A sharper estimate of the data, however, can be obtained using the modified Forchheimer’s equation (13). And so we plot in Figure 14, the friction factor $ f K $ as a function of Reynolds number $R e K $ for

*D*

_{σ}= 0.125, 0.25, 0.375 and for porosities of

*σ*= 0.45 (triangle), 0.6 (square), and 0.7 (circle). Also shown as the dashed curve is the modified Forchheimer’s equation (13). From the figure, we see that Forchheimer’s equation predicts the modified friction factor for small values of the Reynolds number, but at the higher Reynolds numbers, the friction factors begin to deviate from the curve. This suggests that the parameter

*F*should be a function of porosity. Using a least squares fit through our data, we find that

This form is akin to that given by Ref. 35, namely,

where *A*, *B*, *C* are empirically fitted parameters that depend on the microstructure. The modified Forchheimer’s equation (13) with *F* given by (14) is plotted as solid curves for the three values of the porosity. Figure 14(b) shows a blow up region at the higher Reynolds number. Note that the proposed Forchheimer equation with *F* given by (14) well predicts the behavior over the values of the porosity, polydispersity, and Reynolds numbers considered in this study. However, absent a rigorous theory to guide in the fitting process, we make no claim that the fit is unique; our only purpose here is to show that a scaling does exist.

Finally, we show in Figure 15, contours of the velocity magnitude for creeping flow and for the inertial regime (*Re* = 51.9) and for *σ* = 0.7 and *D*_{σ} = 0.375. Note the difference in scales of the magnitude.

## IV. SUMMARY

In this paper, we have examined fluid flow through random periodic packs of polydisperse cylinders. The cylinders are taken to be immobile, and a separation distance is used to avoid blockage in two-dimensions. We solve the two-dimensional unsteady Navier-Stokes equations on a Cartesian grid and use the immersed boundary method to treat internal flow boundaries. The flow solver has been validated for various known problems. Our results indicate that the permeability, when normalized by the length scale *L _{p}*, is a function of porosity and polydispersity. In addition, we find that the permeability can be written as a function of the mean shortest Delaunay edge, and that when properly scaled, the results for the polydisperse case essentially collapses onto the monodisperse case. In addition, we carry out a statistical analysis of the permeability computed from 500 samples and show that (i) rare events can exist, and (ii) the distribution is better characterized using a Weibull distribution rather than a normal distribution. Finally, we carry out simulations in the inertial regime. We find that the Forchheimer equation proposed by Ref. 25, with the parameter

*F*given by Eq. (14), well characterizes the flow over the range of Reynolds numbers considered here. The long term goal of this work is to examine flow through three-dimensional random periodic packs of polydisperse particles and to determine correlation laws for global quantities, such as permeability, that may be used in models at the macroscale.

## Acknowledgments

This work was supported in part by the Defense Threat Reduction Agency, Basic Research Award No. HDTRA1-13-0010 to University of Illinois at Urbana-Champaign and HDTRA1-14-1-0031 to University of Florida; Dr. Suhithi Peiris, Program Manager. T.L.J. was also supported in part by the U.S. Department of Energy, National Nuclear Security Administration, and Advanced Simulation and Computing Program, as a Cooperative Agreement under the Predictive Science Academic Alliance Program, under Contract No. DE-NA0002378.