The Permeability of Boolean Sets of Cylinders

Numerical and analytical results on the permeability of Boolean models of randomlyoriented cylinders with circular cross-section are reported. The present work investigates cylinders of prolate (highly-elongated) and oblate (nearly flat) types. The fluid flows either inside or outside of the cylinders. The Stokes flow is solved using full-fields Fourier-based computations on 3D binarized microstructures. The permeability is given for varying volume fractions of pores. A new upper-bound is derived for the permeability of the Boolean model of oblate cylinders. The behavior of the permeability in the dilute limit is discussed. Résumé — La perméabilité de modèles Booléens de cylindres — Ce travail numérique et analytique porte sur la perméabilité de milieux Booléens de cylindres à section circulaire. On considère les deux cas de cylindres très allongés ou très aplatis. L’écoulement du fluide se fait à l’extérieur ou encore à l’intérieur des cylindres. L’équation de Stokes est résolue par méthode de Fourier sur des microstructures binarisées et les perméabilités calculées. Le comportement de la perméabilité d’un modèle Booléen de cylindre au voisinage de la limite diluée est étudié au moyen de bornes analytiques et les résultats discutés.


INTRODUCTION
In the recent years, several so-called "Fourier-based" algorithms have been introduced to solve the problem of Stokes flow in heterogeneous media [1][2][3].These promising methods, which rely on digitized images of the microstructures and do not require meshing, are especially useful for predicting the permeability of materials based on microstructure models or segmented microtomography images.They have been applied to compute the permeability of textiles [4,5], ceramic foams [6] or anode layers in fuel cells [7], making use of 3D models.
Numerical computations on large-scale 3D microstructures are also useful for investigating theoretical and empirical homogenization estimates, like the classical Carman-Kozeny [8,9] formula, Doi's upper-bound [10,11] based on correlation functions or dilute limit expansions [12].The Boolean model of spheres has been considered in the literature [13] and, as expected, follows the Carman-Kozeny approximation [14], whereas the upper-bound of Doi is valid in the dilute limit only, i.e. for a small volume fraction of obstacles.
The permeability of ordered arrays of cylinders has received considerable interest in experimental, theoretical and numerical works [15][16][17][18][19]. Attention has also been devoted in the literature to various random fibrous media, where fluid flows around cylindrical obstacles.Many of the theoretical homogenization methods used for regular arrays, however, such as the construction of trial fields, cannot be employed for random media.To study random media, numerical methods, such as the Lattice-Boltzmann method, and comparisons with analytical formula for regular arrays have been used.In one of the first such study [20], an analytical estimate, based on numerical predictions, is proposed for disordered fibrous media.The estimate has been found to be accurate for porosities less than about 25% [21].A fiber-web material, with a structure similar to that found in fibrous sheets or paper, has been investigated numerically in [22].Extensive numerical computations have been performed in [23] on models of prolate spheroids with varying aspect ratio.More recently, the efficient Lattice Boltzmann method has been used [24] to compute the permeability fibrous media, where the fluid flows outside of the cylinders.In one such study [25], several models made of random cylinders are investigated.An empirical law is proposed for the permeability, based on the theoretical analysis of the permeability of regular arrays of cylinders [26].The effect of the cylinders curvature is found to be small compared to straight cylinders.Computations have also been carried out with varying cylinders aspect ratios.The combined numerical and experimental study in [5] considers fibrous media.The authors investigate in particular the domain of validity of the relevant analytical estimates.In another experimental investigation, Darcy and Darcy-Forchheimer laws have been identified on packings of non-spherical particles, including prisms and cylinders [27].
Much less work has been devoted to the "reverse" model of fluid flow through fibers, despite some early works [28].Systems of random flat cylinders have been investigated in hydrogeological literature on fractured rocks.The permeability of 3D random models of fracture networks has been investigated numerically and theoretically in various works [29,30].Among important geometrical factors that influence the permeability, authors have highlighted the role of contact regions [31] or that of a wide size distribution in fractures [32].
In this work, a systematic numerical investigation of the Boolean model of cylinders, is carried out with emphasis on effective permeability, fluid velocity distributions, and analytical bounds.The study is organized as follows.The microstructure model is introduced in Section 1.1, the equations to solve for Stokes flow are given in Section 1.2.Fast Fourier Transform (FFT) numerical results are presented in Section 2 with permeability of the various models and numerical maps.Analytical upper-bounds are derived in Section 3 for flows in Boolean models of cylinders in the dilute limit.Concluding remarks are given in Conclusion.

PERMEABILITY OF RANDOM MEDIA 1.Boolean Model of Cylinders
The microstructures considered hereafter consist in Boolean sets of cylinders.Boolean models [33] are a classical type of random models which depend on a primary grain and on a Poisson point process [34].In the present work, the primary grain is defined by a distribution of cylinders, and the homogeneous Poisson point process by a spatially-constant intensity n, the mean number of points per unit volume.A cylinder is attached to each point of the Poisson point process, and the Boolean set is the union of cylinders.In the following, we choose to focus on cylinders with high aspect ratio that significantly depart from sphere packings.For such ideal materials, we expect in particular a strong influence of the microstructure on the permeability.
The cylinders have circular cross-section of radius r and height h but their main axis is randomly-oriented.The direction of the axis follows a uniform distribution on the sphere, resulting in (macroscopically) isotropic Boolean sets.Cylinder intersections are allowed so that the volume fraction of cylinder p varies from 0 to 1 and is given by [33]: where V 0 = pr 2 h is the volume of a cylinder.The Boolean sets accordingly depend on n, h and r, or, equivalently q, h and r.In the following, cylinders with very high or small aspect ratio h/r are considered, with p fixed, as detailed below.All microstructures accordingly depend on two parameters, the cylinders volume fraction q and a characteristic size for the cylinder which is either r or h.

Note however that f c
A ? 0 when r/h ?0 as demonstrated for spheroids [35] and ellipsoids [36].Likewise, the cylinders become infinitesimal lines as r/h ?0 and one expects f c B ? 0. Secondly, let h ( r and consider oblate cylinders.The configuration of interest is that of "porous cylinders" (C) where the fluid flows inside cylinders.In the case h ( r, the cylinders have asymptotically the same shape as sections enclosed by two (randomly-oriented) parallel planes.The percolation threshold f c C again tends to 0 (we refer to [37] for a study on the connectivity of such systems).Accordingly, f is defined in the entire interval [f c C ; 1] % [0; 1].Note that, for r = 1, the "reverse" model where fluids flow outside the cylinders presents no interest as the complementary set does not percolate.
Hereafter, all microstructures for models A, B and C are generated on grids of L 3 voxels with L = 512.The cylinders highest dimension is set to 512 voxels in all cases so that h = L for models A, B and r = L/2 for model C. In configurations A and B, we take r = 3 voxels, to keep a large aspect ratio with a cross-section of the cylinders that resembles disks.By comparison, radii lengths of 6 to 12 lattice units were chosen in media containing 100 3 to 300 3 elements in [21].
Likewise, the value of h for model C cannot be as small as 1 due to discretization effects in the numerical method.We take h = 4 voxels for model C. The aspect ratio is accordingly h/r = 171 (A, B) and r/h = 64 (C).Sections of the microstructures used for models A, B and C are shown in Figure 1.

Stokes Flow
We consider a viscous fluid of velocity u satisfying the Stokes equation in the pores: where p(x) is the pressure at point x and l the fluid viscosity.The fluid velocity is zero (u = 0) along the interface with the solid.Periodic boundary conditions are applied on the computational domain, of the form [38]: here # denotes periodicity and DP is the macroscopic pressure drop, oriented along the first axis e 1 of a Cartesian coordinate system.The permeability k is a scalar for isotropic media.It is given by Darcy's law: The above problem is solved using full-field computations with the FFT-Stokes algorithm [3].The number of iterations used to solve each problem ranges from 15,000 to about 35,000, for an error criterion of g = 10 À6 (we refer to [3] for details on the FFT algorithm).The convergence of g toward 0 is roughly a powerlaw of the number of iterations (Fig. 2).In many instances, however, convergence is either faster or slower during the computations.Each computation took about 2 days on a 16-cores machine with a CPU-parallelized code.

FFT RESULTS: PERMEABILITY
The results are compared to the Carman-Kozeny estimate: where c = 5 is an empirical factor and c is the specific surface area.For Boolean models A, B and C, we have [39]: where S 0 /V 0 is the surface/volume ratio of a cylinder.FFT results for the permeability of models A and B are plotted in Figure 3.The data pertaining to model B (flow outside cylinders) is qualitatively similar to the results of Lattice-Boltzmann computations [21] performed for cylinders with aspect ratio h/r % 29, i.e. 6 times smaller than the one used here.Specifically, the authors distinguish three "concentrated" ( f 0.5), "intermediate" (0.5 f 0.8) and "semi-dilute" (0.8 f ) regimes.The intermediate regime corresponds to an exponential scaling law of k with respect  to the porosity k ~À log f.Our data presents a similar scaling law in the region 0.4 f 0.8 with a slower decrease of k as a function of f in the concentrated regime (f 0.3).The FFT results for model B also compare well to the Carman-Kozeny estimates (Eq.5) and (Eq.6) in a quite wide region 0.3 f 0.9 (solid black curve).As expected, the Carman-Kozeny estimate is not valid in the dilute and concentrated regimes, in agreement with [21].For model B, we also make use of the Jackson-James estimate [20], valid for fibrous media in the dilute region f % 1 (solid green curve).
Similar results hold for models A and C (Fig. 3 and 4).Again, FFT data points are close to the Carman-Kozeny estimates (Eq.5) and (Eq.6).This is especially true for model C, except in two narrow regions near the concentrated and highly-dilute regimes.
Histograms of the fluid velocity component u 1 are represented in Figure 5 for the various models.Remarkably, the shape of the distribution takes different forms with respect to the models.This is especially true of models A vs B. The field maps show a very high concentration of the fluid velocity field around selected points, for all models and all values of the porosity, except in dilute regimes.Maps of the velocity field at selected porosity are shown in Figure 6.

UPPER-BOUND FOR THE PERMEABILITY OF BOOLEAN SETS OF OBLATE CYLINDERS
This section is devoted to upper-bounds for the permeability of a Boolean model of cylinders of the "flat" type, relevant to model C. Recall that flow occurs inside the cylinders in model C. The results in this section are based on exact results previously obtained for the covariogram of cylinders [40].
Microstructures with elongated-type cylinders (h ) r), relevant to models A and B, have much more complex covariance functions and are not treated here.
The Doi [10] upper-bound for permeability reads: where s is the specific surface area of the void-solid interface, and F vv , F sv and F ss are the "volume-volume", "surface-volume" and "surface-surface" correlation functions: In the above, Z is the indicator function for the complementary B c of set B (the obstacles), r is the gradient operator and u is a unit vector (|u| = 1).The quantity F vv is given by the probability that two points separated by a distance t lie in B c .The quantities F sv and F ss can also be interpreted in terms of probability, replacing the interface of B by an interphase of vanishingly small width d.The quantity F sv is proportional to the probability (divided by d) that one point lies in B c and the other in the interphase.The quantity F ss is related to the probability that the two points lie in the interphase.Equation ( 7) is a rigorous bound which provides the exact lowest-order correction to the permeability for dilute concentration of spherical obstacles q ?0 [12].The correlation function F vv of a Boolean model of cylinders is given by [39]: where K(t; h, r) (not to be confounded with k) is the normalized mean cylinder covariogram.The latter is defined as the volume intersected by two cylinders as follows: where |C| is the volume of the cylinder C and C -t is C translated by a vector of length t.The quantity is averaged over orientations uniformly distributed on the sphere and normalized by the cylinder volume, so that k(t = 0; h, r) 1.

Case r = 1
We first consider the case of infinitely flat cylinders with r = 1 [40].Note that h and q are kept finite in this limit.Accordingly, the Poisson intensity h (mean number of cylinders per unit of volume) should tend to 0 in a manner inversely proportional to the cylinder volume.Indeed h is related to q by [39]: and so h $ 1/r 2 ?0. Accordingly, the microstructure should be understood as a set of randomly-oriented, very large cylinders, with a very low density (in number per unit volume), so that the cylinders volume fraction remains finite.This limit microstructure could be usefuly compared to Poisson tessellations (or dilated versions thereof).We refer to [41,42] for theoretical results on the latter.The covariogram reads: The equation for F sv and F ss in Boolean media are similar to that of F vv [43,44] but quite more complex.The formula involve first and second-order derivatives at h 0 = h of the following functional: where, C 0 is a thick plane of width h 0 .Straightforward computations lead to: The three expressions for the correlation functions allow one to compute the integrand in Equation ( 7): and so the integral in Equation ( 7) diverges and the upperbound is infinite.The bound accordingly provides no information on the permeability of materials with "infinitely" flat-cylinders.This negative result will nevertheless prove useful for interpreting the results of Section 3.2.

Case r Finite
We now seek for the computation of upper-bounds on the permeability for cylinders of the oblate type but with finite radius r.The isotropized covariogram K of general cylinders with finite r and h was studied in [40,45].The exact form provided in [40] is a rather complicated analytical expression involving incomplete Elliptic functions.In the present work, we use the following approximate expression [40]: In the above, cos À1 is the inverse cosine (also denoted arccos) function.For r > h, (Eq.18) is a very good approximation of K(t) which becomes asymptotically exact when r/h ? 1.More precisely, the maximal difference = sup t |K(t) À K app (t)| is attained when t ' h for all r > h.For fixed h, the maximal error is a decreasing function of r, and ?0 when r ? 1.As an indication, we have % 1.1% for r = h and % 0.1% for r = 4h (see [40] for details).
The derivation of Doi's bounds for general cylinders requires one to compute the volume of the intersection of two cylinders of different heights and radii averaged over uniformly-distributed directions, which we do not have.
Accordingly, hereafter we use the simpler Berryman-Milton bound [12,46] which reads: For cylinders with r > h, replacing F vv with Equation (11) and using Equation ( 18) yields, after integrating: The upper-bound is only useful in the dilute limit.In fact, the latter largely overestimates the FFT data for the permeability and represented in Figure 4 (not shown on the plot).Its interest lies in the dilute regime q ?0. Fix r and expand the above to zero-order correction O(1) in the dilute limit q ?0 and let r ? 1 afterwards.We obtain: At the lowest-order correction, one recovers an asymptotic regime $ 1/[q(log q) 2 ], a higher correction than that predicted by the James-Jackson model [20] which gives $ À1/[q log q] as first-order correction.The latter is itself higher than the dilute expansion for spheres $1/q.
Interestingly, the second leading-order term in Equation ( 21) scales as 1=ð ffiffi ffi q p log qÞ.A similar non-analytic dependence appears in the dilute limit expansion of the permeability of a bed of totally impenetrable spheres, obtained by self-consistent methods [47]: where a is the sphere radius.The ffiffi ffi q p term is related to hydrodynamic screening effects [48,49].For a dilute packing of spheres, as shown in [12], no screening term is predicted by the bound Equation (19) or by Doi's bound [10].As shown in expansion in Equation ( 21), a screening term is indeed predicted by the variational bound, for the particular highlyelongated shapes considered here.Note that in the dilute limit considered here, the obstacles are not cylinders, but have a variety of shapes with very large surface/volume ratio.
Although the bound in Equation ( 19) is in general less accurate than Doi's bound which uses the surface-surface correlation information [12] both are relevant in the dilute limit.Doi's bound predicts the correct leading-order term in the dilute limit for Boolean spherical obstacles.The bound in Equation (19) predicts the correct scaling law $ a 2 /q, with a slightly overestimated prefactor 4/15 % 0.27 instead of 2/9 % 0.22 for the Boolean model of spheres.
Compute now the limit r ? 1 and afterwards q ?0 in Equation (20): The expansion above depends on r unlike Equation ( 22) and tends to +1 which explains the result obtained previously for Doi's bound with r = 1.When taking the limits r ? 1 and q ?0, different regimes are obtained depending on the order one takes the limits.The two regimes appear when the permeability is plotted as a function of q in log-log plot (Fig. 7).The regime change occurs at some points r % r c (q), q % q c (r) which are obtained by equaling the two highest-order corrections Equation (21) and Equation ( 23): When r ) r c (q), or equivalently q =( q c (r), the permeability follows expansion (Eq.23).Expansion (Eq.22) holds in the domain r ( r c (q) or q =) q c (r).

CONCLUSION
In this work, the permeability of various Boolean models of cylinders with very high aspect ratio have been computed, with fluid flow inside or outside cylinders.The flat and elongated cylinders are tantamount to models of planes or lines of finite width.Good agreement is found with the Carman-Kozeny approximation in the intermediate (non-dilute, non-concentrated) regime.The estimate is particularly relevant for the model of plane or flat cylinders.
The behavior of certain upper-bound in the dilute limit has been computed for the model of flat cylinders in the dilute limit, making use of exact formula for the covariogram of cylinders.The upper-bound predicts two distinct regimes depending on the diameter of the cylinders.In the stronglydilute regime, a scaling law $ 1/[q(log q) 2 ] is predicted, whereas a much weaker regime $ log q arises when r ) h/q.Numerical computations, which are quite challenging with the methods currently available, should be carried out to confirm or infirm the predicted behavior.
Let first r ( h and consider a model of prolate cylinders.Two configurations are possible: "porous cylinders" (model A): the fluid flows inside cylinders; or "cylindrical obstacles" (model B): the fluid flows in the complementary set of the cylinders.Let f be the volume of the phase where the fluid flows, so that f = p for model A and f = 1 À p for model B. Macroscopic flow is permitted if the cylinders percolate (A) or if the complementary set of the cylinders percolates (B).The range of useful values for f therefore reads f c A f 1 and f c B f 1 where f c A , f B are the two percolation thresholds for the Boolean model of cylinders.

Figure 1 2D
Figure 1 2D sections of Boolean models of cylinders (black) of the prolate a) and oblate b) types.a) Both phases percolate.Fluid flow occurs inside the black phase (model A) or inside the white phase (model B). b) The black phase only percolates.Fluid flow occurs inside the black phase (model C).

Figure 2
Figure 2 Error criterion of the MINRES algorithm used in the FFT scheme for Stokes flow computation in model C, as a function of the number of iterations, for various values of the porosity.

Figure 3 Figure 4 Figure 5 Figure 6 2D
Figure 3Effective permeability k vs porosity for a Stokes flow inside and outside a Boolean model of highly-elongated cylinders.Symbols: FFT results, solid lines: Carman-Kozeny estimates.