Permeability of fluid flow through a periodic array of cylinders

Accepted Manuscript Permeability of Fluid Flow through a Periodic Array of Cylinders Kannanut Chamsri, Lynn S. Bennethum PII: DOI: Reference:

S0307-904X(14)00286-8 http://dx.doi.org/10.1016/j.apm.2014.05.024 APM 10031

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

30 May 2013 24 March 2014 16 May 2014

Please cite this article as: K. Chamsri, L.S. Bennethum, Permeability of Fluid Flow through a Periodic Array of Cylinders, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm.2014.05.024

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Permeability of Fluid Flow through a Periodic Array of Cylinders Kannanut Chamsri1 King Mongkut’s Institute of Technology Ladkrabang, Bangkok 10520, Thailand

Lynn S. Bennethum2 University of Colorado Denver, Denver, CO, 80202

Abstract The three-dimensional model of an incompressible Newtonian viscous fluid slowly flowing through a periodic array of cylinders is considered. We use homogenization to determine a system of equations that are then solved numerically to calculate the permeability. We determine the permeability as a function of the angle the cylinders make with the bottom surface. Numerical results are obtained using a Taylor-Hood mixed finite element method. The numerical approach is validated by comparing the results with computed results of Rodrigo P.A. Rocha and Manuel E. Cruz, for a simple cubic array of spheres, with good agreement. Results are presented for different cylindrical densities and angles. When the flow aligns and is perpendicular to an array of cylinders, the results are validated with experimental data. The spherical part of the permeability is compared with the Kozeny-Carman equation. Applications of these results include modeling fluid flow through biological hairlike structures such as animal hair, glass rods, fiberglass, filter pads, polymer gel, collagen, nylon fibers and natural rice field or trees. Keywords: Permeability, Conductivity, Array of Inclined Cylinders, Finite Element Method, Kozeny-Carman equation

Email addresses: [email protected] (Kannanut Chamsri), [email protected] ( Lynn S. Bennethum) 1 Phone number: +66900100827 2 Phone number: (303) 556-4810

Preprint submitted to Applied Mathematical Modelling

May 27, 2014

1. Introduction There are many examples where one may want to analyze the flow of a Newtonian fluid (gas or liquid) through an array of parallel cylinders, including flow through animal hair, glass rods, fiberglass, filter pads, polymer gel, collagen, nylon fibers, grassy vegetation and biological hairlike structures. Modeling these structures as thin porous materials with a solid phase composed of a periodic array of cylinders requires calculating the permeability as a function of the angle that the cylinders make with the bottom surface. Take for example the flow of water through a vegetated wetland. Understanding the drag caused by grassy vegetation has important ramifications in the prediction of water levels for wetland protection and restoration and storm surges on low-lying coastal areas [1]. As water moves, the grassy vegetation bends to varying degrees, requiring the permeability of the grass for a range of angles. The goal of this paper is to do a numerical study to calculate the permeability for incompressible fluid flowing through a periodic array of cylinders as a function of cylinder density and cylinder angle. For an array of parallel cylinders, there are many experimental [2, 3, 4, 5, 6, 7, 8, 9], analytical [10, 11, 12, 13] and numerical [1, 14, 15] results for the permeability or drag for flow both perpendicular and parallel to an array of cylinders. Some of these results pertain to calculating the drag coefficient as oppose to the permeability but the two results are related. While the drag coefficient is used to quantify the resistance acting on an object due to the flow of fluid around it, the permeability measures the resistance fluid has while moving through a porous medium. To the authors knowledge this is the first paper to calculate the permeability tensor for a periodic cylindrical medium where the cylinders form an arbitrary angle to the horizontal plane. To obtain the permeability tensor, we choose to use numerical homogenization, which is generally well-accepted for upscaling a linear equation such as Stokes flow see e.g. [16, 17] and [18]. We note here that upscaling a nonlinear equation (not considered here) requires special consideration and we refer the reader to [18] for an example of upscaling the nonlinear equation representing a non-Darcy flow model. The geometry of our problem is presented in Section 2. In Section 3 we apply the homogenization method to the periodic cell array of cylinders to obtain the system of equations that is discretized in Section 4 using a mixed finite element method. In Section 5 we validate the numerical results by observing their expected behaviors, comparing with numerical results of 2

ΓS

Ω

r

ΓF

r

F

d

lc

Ωs

d

Figure 1: Geometry: an ideal cell of cylinders when the angle between cylinders and horizontal plane is π/2.

Rocha & Cruz [19] for a simple cubic array of spheres, and comparing with experimental data. In Section 5.3, since experimental results from Sullivan [4] were expressed as drag coefficients instead of the permeability, we first explain the conversion of drag coefficient to permeability. Moreover, four different sets of experimental data from three authors [4, 6, 7] are chosen to compare with our numerical values in Section 5.3.2. A polynomial approximation of each entry of the numerical permeability tensor is provided in Section 6. Because the Kozeny-Carman equation is a popular equation to determine the permeability, we compare 1/3 the trace of the permeability of the numerical results with the Kozeny-Carman equation in Section 7. In Section 8 we draw conclusions. 2. Geometry and Assumptions In this work, an array of cylinders is assumed to align periodically along a surface. Figure 1 illustrates a periodic array of cylinders that make a right angle with the bottom surface. The left figure is a three-dimensional ideal cell of cylinders, while the right figure shows the geometry as viewed from the top. The variable r is the radius of the cylinder and d is the distance between the cylinders. The relationship between d and r is d = lc − 2r where lc is the length of the square domain in the right figure. The angle θ between the cylinders and the horizontal plane is illustrated in Figure 2. 3. Homogenization In this section we present the equations that are solved to determine the permeability. Homogenization is a method used to upscale governing equa3

θ Figure 2: the angle θ between the cylinder and the horizontal plane.

tions e.g. Stokes equation, from the microscopic to the resulting macroscopic scale, and has been widely used since the 1970s [20]. One of the strengths of this approach in this study is that for a given microscopic geometry, the coefficients in the resulting macroscopic equations can be found by using periodicity. In particular, by using homogenization to upscale the Stokes equation we get a system of equations that can be used to determine the permeability. In this section, we summarize work from Sanchez-Palencia & Zaoui [21]. We assume that readers have some familiarity with homogenization and only summarize the key assumptions. For more details, see e.g. Bensoussan et al. [20] and Sanchez-Palencia & Zaoui [21]. Let Ω be the domain which consists of a fluid domain, ΩF , and a solid domain, ΩS , with a smooth liquid-solid interface Γ = ΓS ∪ ΓF where ΓS and ΓF are the boundaries of the solid and liquid phases respectively (see Figure 1). We assume slow fluid flow, a rigid solid phase and a viscous incompressible fluid, so that Stokes equation is applicable: ∇·v=0 −∇p + µ 4 v = 0 v=0

in ΩF in ΩF on ΓS ,

(1) (2) (3)

where equation (3) represents a no-slip boundary condition on ΓS ; v denotes the fluid velocity and p and µ are the fluid pressure and dynamic viscosity, respectively. Let L be a macroscopic length scale, typically on the order of at least 100 times lc . We assume the diameter of the cylinders, a = 2r, is small compared to a macroscopic length scale L, i.e. if  = a/L, then   1. Let x be a macroscopic variable and define the microscopic variable or the stretched coordinate y = x/. The dynamic viscosity coefficient µ is assumed fixed and independent of . We consider an asymptotic expansion of v and

4

p, in the form, v = 2 (v0 (x, y) + v1 (x, y) + ...) and p = 0 (p0 (x, y) + p1 (x, y) + ...)

(4) (5)

where vi and pi are Ω-periodic in y, and the exponents on  are chosen to yield a physically meaningful solution. Substituting expansions (4)-(5) into the incompressible system of Stokes equations (1)-(3) and applying the homogenization method, we extract the following system of equations that can be used to determine the permeability: ∇y · ui = 0 i

i

−∇y q + 4y u + e = 0 u=0

in

ΩF

(6)

in on

ΩF ΓS ,

(7) (8)

where ∇y and 4y represent the gradient and Laplacian operators, respectively, with respect to the microscopic scale y and ui and q i are Ω–periodic. The vector ei is the unit vector in the direction of the yi axis, i = 1, 2, 3. To determine k, we consider the weak form of (6)–(8): Find ui ∈ H such that Z ΩF

∇y ui : ∇y h =

Z ΩF

hi dy

(9)

for all h ∈ H, where H = {h : h ∈ (H 1 (ΩF ))3 , h is Ω–periodic, h = 0 on ΓS and ∇y · h = 0} and then the permeability is given by kij =

1 Z ui dy, |Ω| ΩF j

(10)

where uij is the j th component of ui . Remark 1. It is obvious from (9) that the permeability kij is symmetric when h is replaced by uj . 4. Discretization of the Model Problem In this section, we numerically approximate the solution to equations (6)-(8) using a mixed finite element method. This method is employed to obtain an approximate solution to the system of equations using a variational formulation [22]. Let H = {u : u ∈ [H01 (Ω)]3 } 5

(11)

where u is a periodic function and H01 (Ω) = W21 (Ω) is the Sobolev space with zero at a boundary ∂Ω ⊂ Γ where ∂Ω has a positive measure. Let L = {q ∈ L2 (Ω) :

Z

q dy = 0},

(12)

Ω

where q is a scalar periodic function and L2 (Ω) = W20 (Ω) is the Sobolev space. The weak formulation is: Find (ui , q i ) ∈ H × L such that for all Qi ∈ L and for all wi ∈ H −

Z  Ω

Z  Ω



∇y · ui Qi dy = 0, 

−q i ∇y · wi + ∇ui : ∇wi dy =

Z

(13) ei · wi dy.

(14)

Ω

Once the system of equations (6)–(7) is solved, the permeability tensor is obtained from (10). An alternative expression for evaluating kij is 1 Z ∇uj : ∇ui dy kij = |Ω| ΩF

(15)

obtained by replacing h in (9) with uj and using Remark 1. Equation (15) is used as a consistency check with our numerical results kij from (10). In this work, we use the Taylor-Hood isoparametric 10-node tetrahedra, i.e. M = 10 and N = 4 for quadratic and linear functions, respectively, which is generated by the open software Netgen [23], for the finite element discretization in order to ensure the Ladyzhenskaya-Babuska-Brezzi (LBB) consistency condition is satisfied, [24]. An effective reference for threedimensional finite element programming is Kwon & Bang [25] and good references on programming and implementing periodic boundary conditions can be found in Pask et al. [26] and Sukumar & Pask [27]. 5. Validation of Numerical Results We validate the code by verifying the permeability tensor behaves as it should theoretically, comparing the values with other numerically calculated values and by comparing the results with experimental results. 5.1. Expected Behavior We first check that the permeability tensor is independent of the number of cylinders in the reference cell by calculating the permeability tensor for 6

three different periodic cell array of cylinders. The smallest cell consists of 5 cylinders, see Figure 3, with the radius r = a/2 = 0.1, d = 0.3 and the height h = 0.5, and with cell dimension 0.5 × 0.5 × 0.5. The other two cells contain 13 and 25 cylinders with the same radius r and separation distance, d, and with 13-cylinder cell dimension 1×1×1 and 25-cylinder cell dimension 1.5 × 1.5 × 1.5. The three different cells with mesh generation is shown in Figure 3. As is theoretically expected, the values of the permeability tensors calculated from the three different reference cells only differ up to machine precision when a fine enough mesh is used [28].

Figure 3: Three different periodic cells perpendicular to the horizontal plane tilted for easier viewing with a sample mesh generated by Netgen when the radius r = 0.1 and the distance d = 0.3. The height of each cell from left to right is 0.5, 1 and 1.5, respectively.

In Figure 4-6, we present the graphs of the components of the permeability tensor, which the unit is length2 , as a function of r/d (Figure 4) and as a function angle θ (Figures 5-6). Note that when the angle between the cylinders and horizontal plane is π/2, k11 = k22 as can be seen in Figure 4 demonstrating that geometric symmetry yields appropriate symmetry in the permeability tensor. Moreover, for each ratio r/d, the values of k33 are twice that of k11 or k22 implying that the resistance doubles when the cylinders are orientated orthogonal to the flow direction. This result was first justified theoretically by Jackson & James [29] and reduces computational requirements. Figure 5 shows the diagonal entries of the permeability tensor as a function of distinct angles while Figure 6 show off-diagonal component k13 of the tensor. The values for k12 and k23 are within numerical noise (plus or minus 1×106 ) of zero and so for this paper we give values of zero. Because 7

0.03 k11 k22 k33

permeability, length2

0.025

0.02

0.015

0.01

0.005

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

r/d

Figure 4: Permeability of fluid flow through a cell array of cylinders with increasing r/d, θ = π/2. When r/d = 1.2077, the cylinders are touching.

k is symmetric, only k12 , k13 and k23 are presented. 5.2. Comparison with Numerically Calculated Values for Spheres To further validate the code, we compare our results with those of Rocha & Cruz [19], where the geometry is a simple cubic array of spheres and the solutions are obtained numerically. We modified our code so that instead of a cylindrical geometry, we computed the permeability tensor over the periodic cubic cell, Ω, consisting of a single sphere of unitary diameter at the center of Ω and cell length {π/(6s )}1/3 , where s is the volume fraction of the solid sphere. Because the sphere is isotropic, the permeability tensor is a scalar multiple of the identity, kI. Table 1 shows the non-dimensional permeability for solid volume fraction s = 0.216. The variables ks and kRM = 0.03457 denote the permeabilities in this research and by Rocha & Cruz [19], respectively; δk is the relative error of ks with respect to kRM ; #dof is the number of degree of freedom. Note that as the number of degrees of freedom, which in this case is the number of nodes used to calculate the velocities and pressure increases, the error decreases. With approximately 80,000 degrees of freedom, the relative error is less than 0.7% with CPU time about 4.6 hours. Although the error can be decreased to 0.173% more time is required on the order 3.8 days. Different values of solid volume fractions are compared and the (nondimensional) results are presented in Table 2. For 8

−3

4

x 10

k11 k22 k33

3.5

permeability, length2

3 2.5 2 1.5 1 0.5 0 20

30

40

50

60

70

80

90

angle

Figure 5: Diagonal values of the permeability tensor as a function of angle θ for r/d fixed at 1/3.

−4

8

x 10

k13 7

permeability, length2

6 5 4 3 2 1 0 −1 20

30

40

50

60

70

80

90

angle

Figure 6: Off-diagonal values of the permeability tensor k13 as a function of angle θ for r/d fixed at 1/3.

9

Table 1: Non-dimensional permeability through the simple cubic array of spheres with solid volume fraction s = 0.216 and radius = 0.5 unit; ks denotes the permeabilities in this research; δk is the relative error of ks with respect to kRM which is the permeability calculated by Rocha and Cruz (2009).

#dof CPU time (sec) ks δk(%) 11,329 144.31 0.03519 1.793 81,470 16,804.78 0.03480 0.665 615,812 326,804.18 0.03463 0.173

Table 2: Non-dimensional permeability through the simple cubic array of spheres with varying volume fraction of solid s with radius 0.5 unit; ks and kRM denote the permeabilities in this research and by Rocha & Cruz [19], respectively; δk is the relative error of ks with respect to kRM .

s 0.125 0.216 0.343 0.450

ks 0.1037 0.03463 0.01064 0.004419

kRM δk(%) #dof 0.1036 0.096 619,656 0.03457 0.173 615,812 0.01052 1.140 81,660 0.004398 0.477 78,740

s = 0.125, the relative error is smallest, 0.096%, but this case requires a large number of degrees of freedom, 619,656. Although we didn’t increase the number of degrees of freedom for the larger solid volume fractions, we see from Table 1 that the largest relative error is 1%. The compared results are in good agreement for a simple cubic array of spheres. 5.3. Comparison with Experimental Values In some experimental data, the drag force is calculated instead of permeability. The drag force is normally used when the solid phase is viewed as an impediment to flow whereas permeability is used when the fluid is viewed as flowing through the solid phase. Before we compare our numerical permeability to experimental data in Section 5.3.2, we first show the relationship between the drag coefficient and the permeability.

10

5.3.1. Relationship between Drag Coefficient and Permeability To be able to find the relationship, we follow a concept in drag theory that states the total drag force on an array of cylinders in a cell is the pressure drop across the cell [7], i.e. the pressure gradient in Darcy’s Law can be replaced by the drag force. Using this definition and a formula of total drag force of Lamb [30], Chen [7] adapted concepts used in filters. His equation obtained theoretically is: aa 2s ρ ∆p = Cd va2 2 , l π as

(16)

where ∆p is the pressure drop across a cell array of parallel, vertical, cylinders; l is the thickness of the cell; s is the solid volume fraction; Cd is the drag coefficient to be determined; aa is the arithmetically averaged fiber diameter and as is the fiber diameter obtained from surface area measurements; and va = vs /l where l is the porosity andf vs is the superficial velocity of fluid in a cell. As the fibers become more dense, the drag coefficient is modified via the volume fractions. For a periodic array of cylinders, equation (16) becomes ∆p 2s ρ 2 = Cd v (17) l πa a where we assume aa = as = a. Thus aa /a2s = 1/a where a is the diameter of the cylinder. Using va = vs /l as above, we have 2s ρ va ∆p = Cd vs . (18) l πa l Applying Darcy’s Law to (18) and, for simplicity, assuming all variables are scalars, the permeability and drag force coefficient are related by k=

π(1 − s )a2 1 2π(1 − s )r2 1 µπa 1 l = = , 2s ρ Cd va 2s Re Cd s Re Cd

(19)

where Re = ρva a/µ is the Reynolds number and r is the radius of the cylinder. Hence, k∗ =

1 2π(1 − s ) , Cd s Re

(20)

where k ∗ = k/r2 is the dimensionless permeability. Equation (20) will be used to calculate the dimensionless permeability from the drag coefficient Cd . 11

5.3.2. Comparison with Experimental Results We now compare our results with two different sets of experimental data for an array of cylinders: Sullivan [4] where the flow is parallel to an array of cylinders and Brown [6] and Chen [7] where the flow is perpendicular to the array of cylinders. Sullivan [4] presents experimental results for flow parallel to cylindrical fibers such as glass wool, blond and chinese hair, and goat wool. The data from Chen [7] is given in terms of the drag coefficient. For this reference we use (20) to find the permeability. Figure 7 shows the ∗ plots of our dimensionless numerical results k33 and experimental data from Sullivan [4] displaying the harmonization. Similarly, when θ = π/2, the data from Brown [6] and Chen [7] are compared with our dimensionless numerical ∗ ∗ results in Figure 8. The data compare well with ((1/2)(k11 + k22 )). Note that the data from Chen [7] oscillate because he collected data from several publications.

Sullivan (glass wool, blond and chinese hair) permeability

2

Experiment k*33

1.5 1 0.5 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

solid volume fraction Sullivan (goat wool) permeability

50

Experiment k*33

40 30 20 10 0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

solid volume fraction ∗ , of fluid flow parallel to a periodic array of Figure 7: Dimensionless permeability, k33 cylinders for θ = π/2. Figures show the experimental scalar permeabilities of Sullivan [4] for hairs and goat wool, respectively, compared with our numerical results k33 .

12

Brown(Glass Wool) permeability

2

Experiment (1/2)(k*11+k*22)

1.5 1 0.5 0 0.05

0.1

0.15

0.2

0.25

0.3

solid volume fraction Chen (Filter pads) permeability

80

Experiment (1/2)(k*11+k*22)

60 40 20 0

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

solid volume fraction Figure 8: Dimensionless permeability of fluid flow perpendicular to a periodic array of cylinders for θ = π/2. Figures are comparing the experimental data of Brown [6] and ∗ ∗ Chen [7] with our numerical (1/2)(k11 + k22 ) when the flow is perpendicular to the array of cylinders.

13

6. Numerical Permeability Functions In this section, we provide polynomial approximation of each entry of the permeability tensor as a function of the dimensionless distance between the cylinders and the angle the cylinders make with the base. To get an idea of how the permeability varies with the geometry, consider Figures 9. The top left graph presents the numerical permeability component k11 as a function of the dimensionless quotient between the radius of a cylinder, r, and the distance between them, d (x-axis); and the angle θ (y-axis) between cylinders and the horizontal plane defined in Figure 2. The variable r takes on values between 0.025 and 0.125 and the angle θ takes on values between arctan(0.5) to arctan(inf) or about 26 to 90 degrees. Next, we approximate each permeability component by a fourth-order polynomial of the form a1 x4 + a2 x3 y + a3 x3 + a4 x2 y 2 + a5 x2 y + a6 x2 + a7 xy 3 + a8 xy 2 + a9 xy + a10 x + a11 y 4 + a12 y 3 + a13 y 2 + a14 y + a15 , where the coefficients ai , i = 1, 2, ..., 15 are provided in Tables 3. In the right top graph of Figure 9, 4th -order polynomial approximation of k11 is presented. The absolute difference between the calculated values of k11 and the 4th -order polynomial is presented in the bottom graph of Figure 9. Similar graphs can be obtained for the other components. Recall that k12 and k23 are set to zero. The L2-norm for the difference between the numerical results and polynomial approximations for each component is presented in Table 4. 7. Comparison with Kozeny-Carman Equation The Kozeny-Carman equation is an empirical relationship between the porosity and the isotropic, scalar, permeability [31], (l )3 g , k= k1 S 2

(21)

where k = (1/3)tr(k). Equation (21) is derived from the relationship between the superficial or Darcy velocity and the pressure gradient u=

(l )3 g 4P , k1 µS 2 L

(22)

where k1 is a constant depending on the geometry and S is the surface area of the solid phase divided by the total volume of the cell (interfacial area per volume of cell), [S] = [1/L]. Because this expression is so widely used 14

k11

polynomial approximation of k11 0.015 permeability

permeability

0.015 0.01 0.005

0.01 0.005

0 2

0 2 1.5

1.5

0.4 1 0

0.4 1

0.2 0.5

angle

0.2 0.5

angle

r/d

0

r/d

The pointwise error and the twoïnorm = 2.543eï04 ï4

x 10 1.5

Error

1 0.5 0 2 1.5 1 0.5

0

angle

0.1

0.2

0.3

0.4

0.5

r/d

Figure 9: The top left graph is the numerically generated permeability component k11 while the top right one is the fourth-order polynomial approximation of the top left graph; the pointwise absolute error are shown at the bottom.

15

Table 3: The fourth-order polynomial functions: a1 x4 + a2 x3 y + a3 x3 + a4 x2 y 2 + a5 x2 y + a6 x2 + a7 xy 3 + a8 xy 2 + a9 xy + a10 x + a11 y 4 + a12 y 3 + a13 y 2 + a14 y + a15 approximating k11 , k22 , k33 and k13 ; k12 = k23 = 0; x = r/d; y = θ in radius which is defined in Figure 2; r ∈ [0.025, 0.125] and θ ∈ [arctan(0.5), arctan(inf)], about 26 to 90 degrees.

Coefficients a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15

k11 k22 k33 k13 1.0198e + 00 9.8411e − 01 1.3062e + 00 2.9766e − 01 −1.4655e − 02 −1.4485e − 01 −4.6469e − 01 7.0222e − 02 −1.3507e + 00 −1.0963e + 00 −1.2944e + 00 −4.5927e − 01 −6.0732e − 02 −4.2480e − 02 −7.8887e − 02 −1.1961e − 01 1.3877e − 01 2.5436e − 01 6.8503e − 01 1.7072e − 01 6.7157e − 01 4.2544e − 01 3.6564e − 01 1.6761e − 01 −2.8475e − 02 6.2041e − 04 5.2514e − 02 −1.9695e − 03 1.4564e − 01 4.1038e − 02 −9.7006e − 02 1.0366e − 01 −2.1774e − 01 −1.6277e − 01 −1.9823e − 01 −1.7954e − 01 −1.1215e − 01 −5.0873e − 02 −3.4927e − 02 1.1687e − 02 2.7875e − 03 −2.6064e − 05 1.9048e − 03 6.2542e − 03 −1.9444e − 03 −1.1100e − 03 −3.2066e − 02 −2.6424e − 02 −2.6520e − 02 −7.1351e − 03 7.3206e − 02 1.9033e − 02 5.1857e − 02 3.2474e − 02 −1.5704e − 02 1.4396e − 02 8.0453e − 04 −2.7328e − 03 8.2408e − 03 −4.1563e − 03

Table 4: L2-norm errors of the components of the permeability k

components k11 k22 k33 k12 k13 k23

L2-norm error 2.54e − 04 1.32e − 04 2.48e − 04 8.69e − 06 2.11e − 04 6.86e − 07

16

[32, 33, 34], we compare our numerical results with this relationship, even though the isotropy assumption for an array of parallel cylinders is clearly not valid. For the parallel cylinders, there is a one-to-one relationship between the angle and porosity that can be determined analytically. To determine the scalar permeability for the anisotropic case, we substiq 2 l tute the expression S = (1 −  ) 2(sin θ + 1)/(r sin θ) derived in Chamsri [28] into (21) to obtain k=C

(l )3 r2 sin2 θ , 2(1 − l )2 (sin2 θ + 1)

(23)

where C = g/k1 is a constant. Since g is the gravity and the constant k1 depends on the geometry, C also depends on the geometry. Carman [31] determined k1 experimentally when fluid flows through granular beds but not cylindrical ones. To find the constant C for our cylindrical case, we use least-squares regression to get C = 2.0128. Then,

k=

2.0128 (l )3 r2 sin2 θ , 2(1 − l )2 (sin2 θ + 1)

(24)

which is shown in Figure 10 for radius r = 0.1. This is not a bad fit given the anisotropic nature of the cylindrical system. 8. Conclusion In this paper, the permeability of fluid flow through an array of parallel cylinders has been calculated numerically using homogenization and a finiteelement method with a periodic cell domain. This approach indicates that we can find the permeability that is numerically inexpensive. The numerical results were validated by comparing with numerical and experimental results of previous works with good agreement [4, 6, 7, 19]. The diagonal components of the permeability decreases as the angle declines for fixed r/d, see Figure 5. For a constant angle, the diagonal components of permeability decreases with increasing r/d which is consistent with our physical expectation, see Figure 4. In addition, we show that the Kozeny-Carman equation can be a reasonable approximation for 1/3tr(k). Finally, we provide a fourth-order polynomial that can be used to approximate the components of k as a function of r/d and θ, see Table 3. 17

−3

2.5

x 10

our numerical result, tr(k)/3 k(Kozeny−Carman)

permeability

2

1.5

1

0.5

0 0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

porosity Figure 10: The dashed line represents the the spherical part of the permeability tensor of the numerical result while the solid line is from Eq. (24), where x-axis is the porosity varied by changing the angle.

18

References References [1] S. A. Mattis, C. N. Dawson, C. E. Kees, M. W. Farthing, Numerical Modeling of Drag for Flow through Vegetated Domains and Porous Structures, Advances in Water Resources 39 (2012) 44–59. [2] P. C. Carman, The Determination of the Specific Surface of Powders, I, Society of Chemical Industry (Trans. and Communications) 57 (1938) 225–234. [3] E. J. Wiggins, W. B. Campbell, O. Maass, Determination of the Specific Surface of Fibrous Materials, The Canadian Journal of Research 17 Sec. B (1939) 318–324. [4] R. R. Sullivan, Specific Surface Measurements on Compact Bundles of Parallel Fibers, Journal of Applied Physics 13 (1942) 725–730. [5] O. P. Bergelin, G. A. Brown, H. L. Hull, F. W. Sullivan, Heat Transfer and Fluid Friction During Viscous Flow Across Banks of Tubes–III: A Study of Tube Spacing and Tube size, Transactions of the ASME Journal (1950) 881–888. [6] J. C. J. Brown, Determination of the Exposed Specific Surface of Pulp Fibers from Air Permeability Measurements, TAPPI 33 (1950) 130–137. [7] C. Y. Chen, Filtration of Aerosols by Fibrous Media, Chemical Reviews 55 (1955) 595–623. [8] M. L. White, The Permeability of an Acrylamide Polymer Gel, The Journal of Physical Chemistry 64 (1960) 1563–1565. [9] J. A. Wheat, The Air Flow Resistance of Glass Fiber Filter Paper, Canadian Journal of Chemical Engineering 41 (1963) 67–72. [10] J. Happel, Viscous Flow Relative to Arrays of Cylinders, AIChE 5 (1959) 174–177. [11] H. Hasimoto, On the Periodic Fundamental Solutions of the Stokes Equations and Their Application to Viscous Flow Past a Cubic Array of Spheres, Journal of Fluid Mechanics 5 (1959) 317–328. 19

[12] S. Kuwabara, The Forces Experienced by Randomly Distributed Parallel Circular Cylinders or Spheres in a Viscous Flow at Small Reynolds Numbers, Journal of the Physical Society of Japan 14 (1959) 527–532. [13] L. Spielman, S. L. Goren, Model for Predicting Pressure Drop and Filtration Efficiency in Fibrous Media, Environmental Science & Technology 2 (1968) 279–287. [14] F. J. Alcocer, P. Singh, Permeability of Periodic Arrays of Cylinders for Viscoelastic Flows, Physics of Fluids 14 (7) (2002) 2578–2581. [15] T. D. Papathanasiou, The Hydraulic Permeability of Periodic Arrays of Cylinders of Varying Size, Journal of Porous Media 4 (2001) 323–336. [16] P. Renard, G. de Marsily, Calculating equivalent permeability: a review, Advances in Water Resources 20 (5) (1997) 253–278. [17] Y. Chen, L. J. Durlofsky, Adaptive Local-Global Upscaling for General Flow Scenarios in Heterogeneous Formations, Transport in Porous Media 62 (2006) 157–185. [18] C. R. Garibotti, M. Peszy´ nska, Upscaling Non-Darcy Flow, International Journal of Computational Fluid Dynamics 7 (1996) 23–48. [19] R. P. A. Rocha, M. E. Cruz, Calculation of the Permeability and Apparent Permeability of Three-Dimensional Porous Media, Transport in Porous Media. [20] A. Bensoussan, J. L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, North-Holland Publishing Company, 1978. [21] E. Sanchez-Palencia, A. Zaoui, Homogenization Techniques for Composite Media, Springer-Verlag, 1985. [22] J. N. Reddy, D. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, Inc., 1994. [23] J. Sch¨oberl, Netgen, automatic mesh generator (2001). URL http://www.hpfem.jku.at/netgen/

20

[24] H. Elman, D. Silvester, A. Wathen, Finite Elements and Fast Iterative Solvers with applications in incompressible fluid dynamics, Oxford University Press, 2005. [25] Y. W. Kwon, H. Bang, The Finite Element Method using MATLAB, CRC Press LLC, 1997. [26] J. E. Pask, B. M. Klein, P. A. Sterne, C. Y. Fong, Finite-Element Methods in Electronic-Structure, Computer Physics Communications 135 (2001) 1–34. [27] N. Sukumar, J. E. Pask, Classical and Enriched Finite Element Formulations for Bloch-Periodic Boundary Conditions, International Journal for Numerical Methods in Engineering 77 (2009) 1121–1138. [28] K. Chamsri, Modeling the flow of pcl fluid due to the movement of lung cilia, Ph.D. thesis, University of Colorado Denver (2012). [29] G. W. Jackson, D. F. James, The Permeability of Fibrous Porous Media, The Canadian Journal of Chemical Engineering 64 (1986) 364–374. [30] H. Lamb, Hydrodynamics, 6th Edition, Cambridge University Press, London, 1932. [31] P. C. Carman, Fluid Flow through Granular Beds, Transactions of the Institution of Chemical Engineers 15 (1937) 150–166. [32] J. G. Berryman, S. C. Blair, Kozeny-Carman Relations and Image Processing Methods for Estimating Darcy0 s Constant, Journal of Applied Physics 62 (1987) 2221–2228. [33] F. L. Buchholz, Model of Liquid Permeability in Swollen Composites of Superabsorbent Polymer and Fiber, Journal of Applied Polymer Science 102 (2006) 4075–4084. [34] W. D. Comper, O. Zamparo, Hydraulic Conductivity of Polymer Matrices, Biophysical Chemistry 34 (1989) 127–135.

21