Abstract
Viscous flow through fibrous media is characterized macroscopically by the Darcy permeability (K _{D}). The relationship betweenK _{D} and the microscopic structure of the medium has been the subject of experimental and theoretical investigations. Calculations of K _{D} based on the solution of the hydrodynamic flow at fiber scale exist in literature only for twodimensional arrays of parallel fibers. We considered a fiber matrix consisting of a threedimensional periodic array of cylindrical fibers with uniform radius (r) and length connected in a tetrahedral structure. According to recent ultrastructural studies, this array of fibers can represent a model for the glomerular basement membrane (GBM). The Stokes flow through the periodic array was simulated using a Galerkin finite element method. The dimensionless ratio K* =K _{D}/r ^{2} was determined for values of the fractional solid volume (φ) in the range 0.005 ≤ φ ≤ 0.7. We compared our numerical results, summarized by an interpolating formula relating K* to φ, with previous theoretical determinations of hydraulic permeability in fibrous media. We found a good agreement with the CarmanKozeny equation only for φ > 0.4. Among the other theoretical analysis considered, only that of Spielman and Goren (Environ. Sci. Technol. 2: 279–287, 1968) gives satisfactory agreement in the whole range of φ considered. These results can be useful to model combined transport of water and macromolecules through the GBM for the estimation of the radius and length of extracellular protein fibrils.
 fiber matrix
 viscous flow
 computational fluid dynamics
 Darcy permeability
the exchange of water and macromolecules from circulating plasma to the interstitial space is a basic phenomenon that regulates important biological processes (3). Several mathematical models have been developed in the past to describe the hydraulic permeability of the capillary wall, with the aim of identifying the physical forces responsible for this function. The majority of models proposed for the hydraulic permeability of extracellular space within the capillary wall, known as fiber matrix models, are based on the assumption that the extracellular space is filled with fibrous material.
Viscous flow through fibrous media (and through porous media in general) is usually described by Darcy’s law, according to which the average fluid velocity vector U is proportional to the average pressure gradient ∇P
The aim of our study was then to investigate the hydraulic permeability properties of a specific fiber matrix model of basement membrane by numerical simulation of water flow across the fibers and to compare the results with the most popular analytical expressions reported in literature. Extracellular matrix has been usually modeled as a random array of fibers (3). However, recent ultrastructural studies of the basement membrane (GBM) of the glomerular capillary wall show a rather regular fibrillar organization with junctions connecting three or four short fibers (strands) to form a threedimensional meshwork (15, 16,21). We then considered an ordered array of cylindrical strands connected as shown in Fig. 1. We determinedK _{D} by numerical simulation of water flow in a periodic element of the array, under a given pressure gradient. We calculated K _{D} assuming several values of fiber radius and length of the geometrical model and compared the results of the numerical analysis with several analytical expressions proposed in the literature. We also obtained a simple interpolating formula that relates K _{D} to the geometrical parameters of the fiber array. We discuss the applications of our results to the problem of water filtration through the GBM. Since recent studies on glomerular membrane permeability indicate that changes in K _{D}of the GBM are associated with development of renal diseases (7,1113, 29), our model could be used to derive more insight on possible changes in GBM ultrastructure that are responsible for loss of hydraulic permeability. Furthermore, our geometrical model can be applied to other fibrous media, and our results can be used to describe viscous flow in such media.
PROBLEM FORMULATION
Available studies based on high resolution electron microscopy (17, 18,24, 30) suggest that the GBM is composed of a threedimensional meshwork of fibrils with a polygonal structure with pores ranging from 4 to 6 nm in diameter. There are biochemical evidences (1) that collagen type IV macromolecules are flexible rods of about 400 nm in length and that they form the mechanical support of the GBM meshwork. Assembly of collagen IV molecules is due to formation of dimers (with the association of the their COOHterminals) or tetramers (by association of the NH_{2}terminal) (1). There are also evidences that collagen type IV molecules align laterally and twist around each other (33) to form a scaffold into which other components of the GBM are incorporated (such as laminin, heparan sulfate proteoglycans, and others) (31). To mimic the spatial arrangement of fibrils within the GBM, as shown by Moriya et al. (24), we reasoned that the most symmetric junction of four filaments is the tetrahedral structure in which any two adjacent filaments form the same angle (109.5°). We thus considered a tetrahedral periodic arrangement of cylinders of uniform radius (r) and length (L), as shown schematically in Fig. 1, where the periodic unit element consists in the hexagonal prism shown in Fig. 1
A. We calculated the volume occupied by the solid (φ) as
The approach used to calculate K
_{D} was to simulate numerically the flow field inside the described domain under a given average pressure gradient ∇P directed along the zaxis. From the results of the numerical analysis, the average velocity Uwas calculated, and K
_{D} was determined using Eq.1. Numerical simulations for several values of φ, ranging from a very dilute (φ = 0.005) to a very concentrated array (φ = 0.7), were performed. For convenience, results are expressed in terms of the dimensionless permeability
The local Reynolds number (Re) of water moving within the assumed meshwork of fibers was calculated as Re = ρUr/μ, where ρ is the water density, μ is the water viscosity, and U = ‖U‖. For water filtration through the GBM, withr of the order of a few nanometers (17, 18, 24), Re is less than 10^{−7}. In this flow regimen the proportionality between U and ∇P holds strictly; thereforeK
_{D} is a welldefined parameter (independent from Re). We then assumed Re = 0 and considered the stationary Stokes and continuity equations for an incompressible fluid
METHOD OF SOLUTION
A Galerkin finite element solution of the boundary value problem represented by Eqs. 4
and
5
and
710
was obtained using the computational fluid dynamics package FIDAP (version 7.5; Fluid Dynamics International, Evanston, IL), running on an HP 9000/750 workstation. The typical finite element mesh used is shown in Fig. 3. The mesh was built in such a way to make the π_{1} and π_{2} surfaces identical, allowing imposition of the periodic boundary condition (Eq. 9
) on corresponding nodes of π_{1} and π_{2}. The number of nodes ranged from ∼5,000 to ∼7,000 according to different values of φ. We used 8node isoparametric brick elements with linear basis functions for the velocity and with a discontinuous constant pressure approximation. A consistent penalty formulation was chosen, replacing the continuity Eq. 5
by
To verify that the calculated K * did not depend appreciably on element type and on pressure approximation, we performed simulations with 27node isoparametric brick elements, with quadratic basis functions for the velocity, and with three different pressure approximations (discontinuous trilinear pressure, discontinuous linear pressure on a local basis, and discontinuous linear pressure on a global basis). The differences calculated among the results of these simulations were less than 0.05%. We also tested the effect of varying the penalty parameter and noticed that changing ε from 10^{−12} to 10^{−9} did not affect the calculatedK* appreciably.
VALIDATION OF THE METHOD
To check the reliability of the described computational method, we performed a series of numerical simulations of the incompressible viscous flow perpendicular to the twodimensional equilateral triangular array of parallel cylinders shown in Fig.4 and compared the results with theoretical analysis available in the literature. This problem has been studied analytically and numerically by several authors with good agreement among them and with experimental results (see Ref. 21 for review). One of the most updated analysis is that by Sangani and Acrivos (27), who found a general series solution for the Stokes flow and determined the coefficients of the series by a numerical method. For small values of the fractional solid volume φ, their results for the dimensionless permeability (K* = K _{D}/r ^{2}, where r is the cylinder radius) agree very well with an analytical expression found previously by the same authors (28) and subsequently improved by Drummond and Tahir (6). Edwards et al. (9) solved the NavierStokes problem in the same geometry by a Galerkin finite element method for several values of Re confirming, in the case Re = 0, the results of Sangani and Acrivos (27).
In the present work, we assumed the flow to be directed along theyaxis in the unit element delimited in Fig. 4 by the bold line. Since lines S
_{1} and S
_{2} are symmetry boundaries for velocity and pressure, and lines π_{1}, π_{2}, and π_{3} are periodic boundaries for the velocity, we restricted the computational domain to onefourth of the unit element, as shown in Fig. 4 (hatched region). We considered the boundary value problem represented by Eqs. 4, 5, and
7
and the following boundary conditions
A consistent penalty formulation was used with a penalty parameter ε = 10^{−8} (chosen to obtain a negligible compressibility error). We used fournode isoparametric quadrilateral elements with linear basis functions for the velocity and with a discontinuous constant pressure approximation. Typical dimensionless element size ranged from 2 × 10^{−2} to 9 × 10^{−2}. The dimensionless permeability was obtained calculating the average velocity U under a given pressure gradient ∇P. Variations of ε from 10^{−10} to 10^{−6} produced changes inK * of less than 0.3% for all values of φ considered.
The results of our calculations are summarized in Table1, where the values of K * derived from Ref. 27 are also reported for comparison. The latter values of K * were obtained from the dimensionless dragf = F/(μU ) (F being the drag per unit length on a cylinder) by means of the relation K * = π/(φf ), which follows from force balance on the unit element. As shown in Table 1, there is an excellent agreement between values of K * obtained by the present numerical analysis and those derived from Ref. 27. We also tested the dependence of the calculated K * values on the direction of the driving pressure gradient. To this purpose, we performed two simulations (for φ = 0.4 and φ = 0.1) of the flow directed along thexaxis of Fig. 4, using the same procedure. CalculatedK* values differ from those determined for flow along theyaxis by less than 0.1%, indicating that, for Re = 0,K* is not importantly affected by the direction of the driving pressure gradient in respect to the cylinder arrangement.
RESULTS AND DISCUSSION
Numerical results for the dimensionless permeability K * of the threedimensional array, for several values of φ, are reported in Table 2 and in Figs. 5 and6. We obtained a satisfactory interpolation of these data using the following formula, with three freely adjustable parameters
At variance to the twodimensional geometry, to our knowledge there are no other reports in the literature of direct determinations, either analytical or numerical, of the local flow fields in threedimensional arrays of fibers. The available approaches to calculate K * are either semiempirical or based on the extrapolation to the threedimensional geometry of solutions obtained for twodimensional geometries. In the following, we compare our numerical results with the most widely used of these approaches.
According to the semiempirical equation of CarmanKozeny (see Ref. 16for review), the Darcy permeability of a porous medium is
Some authors applied the CarmanKozeny relation to the basement membrane (2, 3, 23, 26), assuming, at variance to our assumption, that the volume of fiber intersections is negligible. The hydraulic radius calculated with this assumption is r
_{h} = ½r(1 − φ)/φ; thus Eq. 18
gives the wellknown expression
Besides the CarmanKozeny relation, other analytical approaches have been proposed to calculate K * for a threedimensional fibrous matrix. Among them, we considered those of Iberall (19), Spielman and Goren (29), and Jackson and James (21) that are more often cited in the literature dealing with permeability of basement membrane. We applied the expressions proposed by these authors to our geometrical model and compared calculated K
_{D} with our numerical results (see Fig. 6). At variance with the CarmanKozeny relation, these approaches are expected to give good results only for low values of φ, since they are based on two simplifying assumptions that apply only to dilute fibrous media, as follows: 1) that the total drag exerted by a fluid is given by the sum of the drag on individual cylinders, neglecting their interactions; and 2) that the drag on a single cylinder is the linear combination of the drag exerted by the flow parallel and perpendicular to the fiber. On the basis of these assumptions, the dimensionless permeability K * for our model, with an arbitrary orientation of the average velocity U, is given by
Iberall (3, 19) proposed to calculateK*_{t} andK*_{n} by means of relations derived previously by Emersleben (11) and Lamb (22). In this way, assuming Re = 0, one obtains the socalled drag theory relation
Spielman and Goren (29) determined the drag exerted on a cylinder for both parallel and perpendicular flow with an approach based on the Brinkmann equation and obtained implicit equations relatingK*_{t }andK*_{n} to φ, which can be combined with Eq. 22
to give
The third analytical approach we have considered is that of Jackson and James (20), who determined K * for a cubic array of cylinders. As the fraction of perpendicular and parallel fibers in a cubic array and in our model is the same, these authors actually employed Eq. 22
to derive K*, and they used forK*_{t} andK*_{n} the dimensionless permeabilities of a twodimensional square array of parallel cylinders calculated by Drummond and Tahir (6). We adapted this procedure to our tetrahedral array, considering K*_{t}and K*_{n} of a twodimensional equilateral triangular array, also taken from Drummond and Tahir (6). Truncating the series in Ref. 6 to the second order in φ we get^{2}
On the basis of the above comparisons, we can conclude that none of the analytical expressions we considered predicts K * with very good accuracy in the whole range of φ adopted in the present study. The best agreement is shown by Eq. 24 but is still unsatisfactory. Therefore, we propose to use our interpolating formula (Eq. 17 ) to express K * as a function of φ for an ordered array of fibers such as that considered in our model. WhetherEq. 17 applies also to different spatial arrangements of fibers cannot be established on the basis of the present results and could eventually be decided with more complex sensitivity studies. In the case of basement membranes, since the actual ultrastructure is not known in detail, such studies would not be justified at the present moment. Our model is consistent with available limited observations of the GBM ultrastructure. Since the various analytical expressions considered do not well reproduce numerical data for such a model, we suggest to use Eq. 17, together with the definition ofK * given in Eq. 3, to relate Darcy permeability of the GBM to the radius of the extracellular protein fibrils and the fractional solid volume.
As stated previously, experimental estimation of the dependence ofK _{D} on r and φ is technically difficult. There are only few reports on the dimensions of the protein fibers that compose the GBM. The available theoretical and experimental observations indicate that the radius of these fibers (r) should range from 0.8 to 1.0 nm (26). Much higher values have been obtained (r > 6 nm) using electron microscopy techniques (18, 24), but there are reasons to suspect that tissue processing for scanning electron microscopy may influence the dimensions of the fibrillar material observed within the GBM. Similarly direct measurements and theoretical predictions would indicate that fractional solid volume (φ) of the GBM is about 0.1 (26). Using our interpolating formula, we calculated, assuming a value of φ = 0.1, that K _{D} = 1.36 and 2.13 for r = 0.8 and 1.0 nm. These values are close to the values of the Darcy permeability of isolated rat GBM derived by Drumond and Deen (7), on the basis of the experimental results of Daniels et al. (5) and of Robinson and Walton (26), who reported values of K _{D} = 2.7 and 1.8 nm^{2}, respectively. Although there is great uncertainty on these parameters, it appears that our theoretical approach closely reflects independent estimations of the geometrical parameters of the GBM meshwork and its Darcy permeability.
The present analysis can be useful in estimating the dimensions of fibrous material that composes the GBM. Recent in vitro and in vivo studies (4, 5, 13, 26) reported estimates of K _{D} for the GBM, and there are evidences that changes in GBM hydraulic permeability develop in experimental and human glomerular injury (8,10, 12, 13, 15, 32). Thus, it would be interesting to know whether these changes in permeability are related to changes in the spatial organization of the GBM fibrillar proteins or to dimensional changes of the fibrous or both. Since K _{D} depends on two parameters, however, estimates of K _{D} alone do not allow one to calculate both r and φ, but another independent relation is needed. A way to overcome this limitation could be the simultaneous consideration of the transport of water and macromolecules across the GBM. The permeability of a fibrous membrane to a molecule of given size depends on the geometrical parameters of the fibers. In detail, according to Ogston et al. (25), the permeability coefficients of fibrous membranes to neutral macromolecules can be expressed as a function of r and φ (3). Combining this relationship with our expression of the Darcy permeability given by Eqs. 3 and 17 , one can uniquely determine the values of φ and rfrom direct measurements of hydraulic membrane permeability and filtration rates of neutral macromolecules.
Acknowledgments
We thank Nicola Zatelli and Ehab I. Mohamed for help in the preparation of the manuscript.
Footnotes

Address for reprint requests: A. Remuzzi, Biomedical Engineering Laboratory, Mario Negri Institute for Pharmacological Research, Via Gavazzeni 11, 24125 Bergamo, Italy.

Present address of M. Palassini, Scuola Normale Superiore, Piazza dei Cavalieri 7, 56100 Pisa, Italy.

↵1 The FIDAP code does not allow imposition of periodic boundary conditions on the single components of the velocity, but only on the three components simultaneously. Thus, Eqs. 9 and 10 apparently represent four boundary conditions on the inflow and outflow surfaces. Because of the continuity equation, however, only two of the boundary conditions of Eq. 9 are independent. Therefore, three independent conditions are actually imposed.

↵2 The expression forK*_{n} given in Ref. 6 is valid only for values of φ not higher than ∼0.5. In Eqs. 25 and 26 the neglected higher order terms would give a contribution smaller than 1% for φ ≤ 0.5. The same is true for Eq. 27.Thus, Eq. 27 can be safely used for comparison of our numerical results with the analysis of Drummond and Tahir (6) in the whole range of applicability of the latter.
 Copyright © 1998 the American Physiological Society