- Email: [email protected]

Computer simulation of ﬂow through a lattice ﬂow-cell model A.R. Mazaheri

b

a,b

, B. Zerai c, G. Ahmadi a,b,*, J.R. Kadambi a,c, B.Z. Saylor c, M. Oliver a,c, G.S. Bromhal a, D.H. Smith a,d

a National Energy Technology Laboratory, U.S. Department of Energy, Morgantown, WV 26507-0880, United States Department of Mechanical and Aeronautical Engineering, Clarkson University, Potsdam, NY 13699-5725, United States c Case Western Reserve University, Cleveland, OH 44106, United States d Department of Physics, West Virginia University, Morgantown, WV 26506, United States

Received 24 February 2004; received in revised form 14 October 2004; accepted 18 October 2004

Abstract For single-phase ﬂow through a network model of a porous medium, we report (1) solutions of the Navier–Stokes equation for the ﬂow, (2) micro-particle imaging velocimetry (PIV) measurements of local ﬂow velocity vectors in the ‘‘pores throats’’ and ‘‘pore bodies,’’ and (3) comparisons of the computed and measured velocity vectors. A ‘‘two-dimensional’’ network of cylindrical pores and parallelepiped connecting throats was constructed and used for the measurements. All pore bodies had the same dimensions, but three-diﬀerent (square cross-section) pore-throat sizes were randomly distributed throughout the network. An unstructured computational grid for ﬂow through an identical network was developed and used to compute the local pressure gradients and ﬂow vectors for several diﬀerent (macroscopic) ﬂow rates. Numerical solution results were compared with the experimental data, and good agreement was found. Cross-over from Darcy ﬂow to inertial ﬂow was observed in the computational results, and the permeability and inertia coeﬃcients of the network were estimated. The development of inertial ﬂow was seen as a ‘‘two-step’’ process: (1) recirculation zones appeared in more and more pore bodies as the ﬂow rate was increased, and (2) the strengths of individual recirculation zones increased with ﬂow rate. Because each pore-throat and pore-body dimension is known, in this approach an experimental (and/or computed) local Reynolds number is known for every location in the porous medium at which the velocity has been measured (and/or computed). 2004 Elsevier Ltd. All rights reserved. Keywords: Flow cell; Pore-scale model; Muliphase ﬂows; Porous media

1. Introduction The long-term motivation of this study is to provide a better understanding of the fundamentals of single and multiphase ﬂow through porous media. This topic has many environmental and industrial applications including enhanced oil recovery, DNAPL (dense non-aqueous phase liquid) remediation, and CO2 geologic sequestra* Corresponding author. Address: Department of Mechanical and Aeronautical Engineering, Clarkson University, Potsdam, NY 136995725, USA. Tel.: +1 315 268 2322; fax: +1 315 268 6438. E-mail address: [email protected] (G. Ahmadi).

0309-1708/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2004.10.016

tion. The conventional approach is to use DarcyÕs law, and saturation-dependent relative permeabilities for two-phase ﬂows [1,5,8,31]. Studies of single-phase ﬂow through permeable media usually begin with DarcyÕs phenomenological equation [6] of about 150 years ago. According to DarcyÕs law, the total ﬂow through such a medium is proportional to both the pressure gradient along the macroscopic-average ﬂow direction and the cross-sectional area of the medium. Details of the microscopic geometry or of the ﬂow inside the medium, however, are not addressed. At higher ﬂow velocities, the pressure gradient and ﬂow velocity may no longer be proportional. For

1268

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

such cases Forchheimer [14] suggested that DarcyÕs law be modiﬁed to include a second-order term in the velocity. Later, the idea of non-integer exponents for the relationship between pressure gradient and ﬂow rate was introduced by White [36] and generalized by Missbach [24] ; but no agreement on the exact value of the exponent has been found [31]. Hence these attempts to return from the two terms of the Forcheimer equation to a single term have been of limited value. As described by Scheidegger [31], many attempts have been made to correlate DarcyÕs law with the Reynolds number and a friction factor. Virtually all such correlations have foundered on the problem of attempting to deﬁne the Reynolds number in terms of pore diameter for porous media (whether natural or synthetic) that have no simple tubular pores. Others, exempliﬁed by the work of Kozeny [19], Carman [2], and Ergun [9], have turned to theoretical descriptions of laminar and turbulent ﬂow to elucidate the Forchheimer equation [31]. Lattice Boltzmann (LB) approaches currently oﬀer perhaps the most popular and promising theoretical method for single-phase ﬂow through permeable media [3,15,22,38]. An extensive review of LB modeling was provided by Chen and Doolen [3]. Lattice Boltzmann studies of single-phase ﬂow through media of packedspheres or media with randomly placed identical obstacles were reported by Succi et al. [34], Maier et al. [21], Koponen et al. [18], and Pan et al. [25]. Using an extension of work by Fenwick and Blunt [10], Pan et al. [25] developed a single-phase pore-network model with cubic pore bodies connected by square cross-section channels for incompressible single-phase ﬂow of viscous ﬂuids. Although the LB method is a powerful technique, it is computationally demanding. Pan et al. [26] developed a new approach for modeling ﬂow through porous media using the LB method with higher-performance, particularly, in ﬂoating point operations and storage capabilities. As opposed to the molecular-level approach of the LB method, there is a very large body of literature that assumes Poiseuille type ﬂow through the throats of a network of ‘‘pore bodies’’ and connecting, tubular ‘‘pore throats’’ [4,11,12,16,17,20,28,29,32,33]. Tsakiroglou [35] presented an analytical approach for the prediction of non-Darcian creeping ﬂow of shear-thinning (non-Newtonian) ﬂuids in porous media. In this work, the analytical models were then integrated into numerical simulations of single-phase ﬂow through a two-dimensional, square-lattice, pore-network model. However, almost all of this literature treats problems of ﬂow of two (or more) ‘‘phases,’’ that is, either two ﬂuids that are miscible and thus have no capillary pressure, but have ‘‘not yet’’ mixed; or else two ﬂuids that truly are two phases at chemical equilibrium with a well-

deﬁned capillary pressure that is additive with the Poiseuille viscous pressure. It has been shown that such models can successfully produce the special two-phase ﬂow cases of DLA (diﬀusion-limited aggregation) and invasion percolation (IP) [12], and the crossover from capillary ﬁngering to compact invasion [13]. Such models have seldom (if ever) been used for ﬂow of a singleﬂuid, but are capable of treating single-ﬂuid ﬂow as a ‘‘trivial’’ case. Detailed measurements for the velocity ﬁeld have been of considerable interest for understanding of the complex features of ﬂuid ﬂows and for model veriﬁcation. Particle image velocimetry (PIV) is now a wellestablished technique that provides the velocity vector ﬁeld in a plane of ﬂow. In this approach laser pulses illuminate a plane in the test section. The ﬂow is seeded with tracer particles that reﬂect laser light, and the particle positions at each laser pulse are recorded and tracked for evaluating the instantaneous velocity ﬁeld [27]. In recent years, Santiago et al. [30], Meinhart et al. [23] and Devasenathipathy et al. [7] have extended the approach to so-called micro-PIV for making velocity measurements in small devices. A more systematic approach to analysis of ﬂow through porous media is now made possible by a triad of recent developments: (1) network models, (2) computers and computational methods for ﬁnding non-analytic solutions of the Navier–Stokes equation, and (3) microPIV measurements. Brieﬂy: (1) network models replace complex permeable media with geometries that are computationally and experimentally tractable; (2) detailed, pore-level spatial-distributions of pressure gradients and velocity vectors can now be computed for these networks; and (3) these velocity vectors can be measured in the laboratory by PIV techniques for ﬂow through networks whose geometries are essentially identical to those used in the Navier–Stokes computations. This approach permits detailed comparisons of computed and measured velocity vectors at diﬀerent locations and ﬂow rates. For the lattice Boltzmann approach, measurements of velocity vectors and/or solutions of the Navier– Stokes equation oﬀer an experimental check on the computations and/or an opportunity to compare LB results with results from classical ﬂow equations. For two-phase network models that assume Poiseuille ﬂow, the computational and experimental approach of this paper suggest the possibility of ‘‘going back’’ to explore the limitations of the Poiseuille ﬂow assumption for single-phase ﬂow; such explorations could illuminate the limitations of the assumption for two-phase ﬂow, and/or lead to extensions of this paperÕs approach to include two-phase ﬂow through networks. This paper describes a complementary combination of computations and laboratory experiments on singlephase ﬂows in a lattice network. Viscous ﬂow in a typi-

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

1269

cal network lattice of 20 · 20 pore bodies was studied using the Navier–Stokes equation. A three-dimensional unstructured computational grid was generated. Numerical simulations were performed for a range of ﬂow rates, and the details of the velocity and pressure ﬁelds were evaluated. The eﬀective permeability of the cell was estimated and the cross-over from Darcy ﬂow to inertial ﬂow was examined. An identical experimental ﬂow-cell was also fabricated and the velocity vectors in the mid-plane of several pore bodies and their adjacent throats were measured using a micro-PIV technique. The numerical simulation velocity vectors were compared with the measured vectors in the laboratory ﬂow-cell, and generally good agreement was found.

2. Flow-cell description The experimental 64 · 64 mm ﬂow-cell was made of acrylic plastic and consisted of 20 · 20 pore bodies of 2.5 mm diameter and height 1.0 mm, connected on a diamond lattice by 5.0 mm long, square cross-section throats of sides 0.2, 0.6, or 1.0 mm. Each pore body was connected to four connecting throats. The size of each throat was selected randomly from the three sizes. The ﬂow entered the cell through a 2.5 mm tube into a small plenum that connected to the pore bodies. At the exit area, ﬂow from the cell discharged into a small plenum and left the ﬂow-cell through a 2.5 mm outlet tube. Fig. 1 shows a picture of the experimental ﬂow-cell used in the study. Experimental and computational results are reported for pore bodies B, C, and D and the ﬂow throats connected directly to each of these pore bodies. Each of these pore bodies was connected to four

Fig. 2. Schematic geometry of the ﬂow-cell model used in the simulation.

throats of identical width (B: 1.0 mm; C: 0.6 mm; D: 0.2 mm.) A computer model was also developed with geometry identical to that of the experimental cell. A section of the computer model is shown in Fig. 2. The cell design was suggested by ﬂow cells used by Ferer et al. [11–13] and Smith et al. [33], who performed ﬂow simulations of a much diﬀerent type.

3. Numerical simulations A three-dimensional computational grid for simulation of ﬂuid ﬂow through the ﬂow-cell was developed. The pore bodies and the throats were modeled using an unstructured grid of about 400,000 cells. In the process of creating mesh, the model was examined to make sure that at least ﬁve grid-cells were generated inside the smallest width passages. Grid sensitivity was check by increasing the number of grid cells to 700,000, until no noticeable changes in results were observed. Fig. 3 shows the grid schematic of a portion of the threedimensional model of the ﬂow-cell. The Navier–Stokes and continuity equations governing the motion of an incompressible viscous ﬂuid are given as oui oui oP o2 ui q þ uj þl þ qgi ; ð1Þ ¼ oxi ot oxj oxj oxj

Fig. 1. Picture of the ﬂow-cell and the locations for which experimental data were obtained.

ouj ¼ 0. oxj

ð2Þ

1270

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

Fig. 3. Unstructured grid schematic of a portion of the threedimensional ﬂow-cell model.

Here, q is the density of the ﬂuid, ui is the ﬂuid velocity, P is the pressure, and l is the ﬂuid viscosity. The FLUENTTM code (version 6.02) was used, and the ﬂow properties in the cell for a range of parameters were evaluated.

4. Baseline computational results The baseline computational results for a constant injection rate of 11.13 · 108 m3/s of sodium iodide solution are presented in this section. This ﬂow rate corresponds to an injection velocity of 0.025 m/s, which corresponds to a Darcy velocity of about 0.0017 m/s. The sodium iodide solution had a density of 1700 kg/ m3, and a viscosity of 0.0054 kg/m s. These values were used in the computer simulation model. In the computational analyses, the no-slip boundary condition was imposed. The PIV measurements were conducted for pore bodies B, C, and D shown in Fig. 1. These pores were the same ones for which detailed numerical results are provided in Figs. 5–7. Fig. 4 shows the contour plot for variations of the static pressure in a plane at the midsection of the lattice-cell model. Fig. 4 shows a continuous pressure gradient throughout the cell with the total pressure diﬀerence of about 600 Pa across the lattice-cell. The corresponding permeability of the cell is estimated to be about 109 m2. A plot of computed velocity vectors in a plane at the mid-section of pore body B (with four 1 mm equal-size throat connections) is shown in Fig. 5a. This ﬁgure shows that the ﬂow is entering the pore body from the two left-hand throats, which are on injection side of the cell. The maximum speeds in the upper and lower throats around the inlet are, respectively, about

Fig. 4. Static pressure contours at the mid-section of the ﬂow-cell model.

40 mm/s and 15 mm/s. However, when the ﬂuid leaves the pore body, through the two right-hand throats, the ﬂuid velocities become comparable, with the lower throat having a slightly higher velocity than the upper throat, about 30 mm/s and 20 mm/s, respectively. Since the velocities in diﬀerent throats are diﬀerent, several diﬀerent Reynolds numbers may be deﬁned. We used the hydraulic diameter of the throat with maximum velocity and the corresponding average velocity in that throat to deﬁne a characteristic Reynolds number for the pore body. That is Re ¼

V dh m

ð3Þ

Here dh is the hydraulic diameter of the throat with maximum velocity, V is the average velocity in the throat and m is the kinematic viscosity of the ﬂuid. Here for the pore body B shown in Fig. 5, Re 12.5. It should also be noted here that due to the very short length of the throats, the ﬂows in the throats are not necessarily fully developed with a well established parabolic proﬁle. Fig. 5b shows the pressure contour plot across the pore body in the same plane as Fig. 5a. This ﬁgure shows that the pressure diﬀerence is about 10 Pa between inlets and outlets of the pore body. Fig. 5d shows the variation of the velocity magnitude contour in the lower right throat. It is seen that the ﬂow is developing as it enters from the port body into the throat. The peak velocity is in the centerline of the throat and the velocity decreases toward the walls to meet the not slip boundary condition at the wall. Fig. 6a shows the computed velocity vectors plotted at the mid-section of pore body C, with four 0.6 mm equal-size throat connections. This ﬁgure shows that the ﬂow is entering the pore body from the lower lefthand throat with maximum speed of about 35 mm/s and is leaving the pore body from the two top throats

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

1271

Fig. 5. Computed and measured ﬂow conditions at the mid-section of pore body (B) with four 1 mm equal-size throat connections. (a) Computed velocity vector plot, (b) static pressure contours, (c) measured velocity vector plot from micro-PIV and (d) velocity contours in a connecting tube.

with top speeds of 10 mm/s and 25 mm/s. It is also seen that the ﬂow rate is negligibly small in the lower righthand throat. Using the hydraulic diameter of the square cross-section throat, the characteristic Reynolds number in the pore body C is about 6.6. The pressure contour plot across this pore body, shown in Fig. 6b, indicates that there is a small pressure diﬀerence between the pore body and the lower-right throat with negligible ﬂow rate. This ﬁgure also shows a maximum pressure diﬀerence of about 20 Pa between the left inlet and outlet throats. Fig. 7a shows the computed velocity vector plot at the mid-section of pore body D, with four 0.2 mm equal-size throat connections. This ﬁgure shows a very low ﬂuid velocity inside the pore body, while the liquid velocity is about 10 mm/s in all the connecting throats. This ﬁgure also shows that the ﬂow enters the pore body from the two left channels and leaves through the two right throats. In this pore body, the characteristic Reynolds number is about 1.3. Comparing the charac-

teristic Reynolds numbers in pore bodies, B, C, and D, it is found that the characteristic Reynolds number decreased as the hydraulic diameter of the throats decreased. A pressure contour plot across pore body D is shown in Fig. 7b. It is seen that the pressure inside the pore body is almost constant throughout the pore, and the pressure diﬀerence between the inlet and outlet is 50 Pa. Comparison of Figs. 5–7 shows that the pressure drops across the pore bodies increase as the throat sizes decrease.

5. Experimental measurements Particle image velocimetry, PIV, is a technique for measuring ﬂuid ﬂow velocities. Rapidly pulsed planar laser light sheet illuminates the test section, and images of the tracer particles lying in the light sheet are recorded by the CCD camera though the receiving optics. In this technique, the displacements of the

1272

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

Fig. 6. Computed and measured ﬂow conditions at the mid-section of pore body (C) with four 0.6 mm equal-size throat connections. (a) Computed velocity vector plot, (b) static pressure contours and (c) measured velocity vector plot from micro-PIV.

particle images are measured in the plane of the image and used to determine the velocity of the particles in the ﬂow. Micro-particle image velocimetry (PIV) was used to measure the velocity ﬁelds in the throats and pores of the experimental lattice-cell shown in Fig. 1. The PIV equipment was comprised of a 10 bit Kodak ES1.0/ 10bit (1018 · 1008 pixel) CCD camera, an EPIX PIXCI-D frame grabber board, a microscope objective lens, and a 120 mJ/pulse Nd:YAG laser (New Wave Research). The system parameters are listed in Table 1. Two sets of ﬁeld of view (FOV), 2500 lm · 2500 lm and 5000 lm · 5000 lm, were used to measure the ﬂow ﬁelds in the lattice-cell. Measurements were made in the 2.5 mm diameter tube at the test cell exit, and in the 1 mm, 0.6 mm and 0.2 mm wide throats of the B, C, and D pore bodies shown in Fig. 1. Sodium iodide solution was injected into the lattice-cell at various constant rates. Refractive index matching between the test cell material (acrylic) and sodium iodide solution was used to eliminate the laser-light path distortion in the test cell and scattering of light from the throat and pore surfaces. The test-cell exit tube was also enclosed in a square

acrylic block to eliminate light-path distortions. The measured axial velocity proﬁle in the 2.5 mm diameter test cell exit tube was laminar. The ﬂow rate obtained from the measured velocity proﬁle was within ±1.6% of the ﬂow rate obtained by collecting and weighing the produced ﬂuid. The expected uncertainty levels in the velocity measurements are about ±1.8%. The uncertainty associated with the pore-throat scale size measurement is 1%. Additional details of the experimental method may be found elsewhere [37]. The PIV measurements were conducted for pore bodies B, C, and D shown in Fig. 1. These pores are the same ones for which detailed numerical results were provided in Figs. 5–7.

6. Comparison of the PIV experimental data with computer simulations In this section, the numerical simulation results are compared with experimental data obtained by the micro-PIV technique. In the experimental study, a constant ﬂow rate of 11.13 · 108 m3/s of sodium iodide solution, which is identical to the conditions used in the computer

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

1273

Fig. 7. Computed and measured ﬂow conditions at the mid-section of pore body (D) with four 0.2 mm equal-size throat connections. (a) Computed velocity vector plot, (b) static pressure contours and (c) measured velocity vector plot from micro-PIV (vectors corresponds to ﬂow direction).

Table 1 Laser, camera, ﬂuid and seed-particle parameters used in the micro-PIV measurements Laser energy

Wavelength

Maximum frequency

Camera

Laser and camera 120 mJ/Pulse

532 nm

15 Hz

Kodak ES 1.0 CCD

Fluid properties Liquid

Density

Viscosity

Refractive Index

Flow rate

NaI (74%), Glycerin (21%), H2O (5%)

1.7 g/cm3

5.4 CP

1.485

11.13 · 108 m3/s

Seed particle properties Material

Density

Refractive index

Settling velocity

Mean radius size

Silicon carbide

3.65 g/cm3

2.65

3.2 lm/s

2 lm

simulation, was used. The experimental results for diﬀerent cases are plotted in Figs. 5c, 6c, and 7c for comparison with the computational results shown in Figs. 5a, 6a, and 7c, respectively. Fig. 5c shows the velocity vector plot obtained by the micro-PIV technique at the mid-section of pore body B (with four 1 mm throat connections). This ﬁgure shows that the general pattern of the measured velocity ﬁeld is quite similar to the computed one shown in Fig. 5a. In particular, the direction and magnitude of the ﬂow ﬁeld in the throats as predicted by the computer simulation

are in good agreement with the experimental data. The average measured velocity of about 35 mm/s in the upper left inlet is comparable with the 40 mm/s that is predicted in the numerical simulation. Fig. 5c also shows that the measured liquid velocity leaving the right lower throat of pore body B is about 31 mm/s compared to 30 mm/s estimated from the numerical model shown in Fig. 5a. The velocity vector ﬁeld measured by the micro-PIV technique at the mid-section of pore body C with four 0.6 mm connecting throats is shown in Fig. 6c. This

1274

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

ﬁgure shows a negligibly small measured velocity in the lower right throat of pore body C, which is in agreement with the numerical simulation result shown in Fig. 6a. Both the experimental velocity shown in Fig. 6c and the numerical one shown in Fig. 6a show a maximum inﬂow velocity of about 35 mm/s in the lower left inlet throat. The outﬂow velocity of about 11 mm/s measured in the experimentation in the upper left outlet throat is also consistent with the numerical simulation results. Fig. 7c shows the measured velocity vector distribution at the mid-section of pore body D, with four connecting 0.2 mm throats. Here equal-size vectors are shown for the direction of the velocity. This is because the velocity inside the pore body D is very small and a bubble was also present during the experimentation that distorted the velocity ﬁeld inside the pore body. Nevertheless, the general features of the measured ﬂow pattern are comparable with the computed one. In particular, as in the computer simulation results, the measured velocity ﬁeld shows that the ﬂow enters the pore body from the two left throats and leaves through the right two connecting throats.

7. Computational parametric studies In the previous section, the computational model predictions were compared with the micro-PIV experimental data for selected pore bodies in the ﬂow-cell for an inlet velocity of 0.025 m/s and reasonable agreement was found. In this section, the computational model was used to perform parametric studies with a range of ﬂow rates. The ﬂow ﬁelds in the cell were computed for inlet velocities of 0.025, 0.075, 0.125, 0.25, 0.375, and 0.5 m/s. For pore body E, with one 0.2 mm, one 0.6 mm, and two 1 mm throats, Fig. 8 illustrated the computed velocity vector ﬁelds. Here the results for inlet velocities of

0.025 m/s and 0.5 m/s are compared, so that the ﬂow rate for the high speed case is 20 times the baseline ﬂow rate. In Fig. 8, the velocity is scaled relative to its peak value at each ﬂow rate, so that a comparative study can be performed. This ﬁgure shows that the general features of the normalized ﬂow ﬁelds appear quite similar for the two ﬂow rates. The ﬂow enters the pore body from the lower 0.6 mm and 0.2 mm throats and leaves through the top two 1 mm connecting throats. At the higher injection rate, Fig. 8 shows, the high-speed jet from the 0.6 mm throat penetrates more deeply into the body of the cell compared to the low-speed case. In the high ﬂow-rate case a small recirculation region forms on the right-hand side of the pore body, which was absent in the baseline case at 0.025 m/s. As was noted before, diﬀerent Reynolds number may be deﬁned. Here we used Eq. (3) to deﬁne Reynolds number for the pore body. For injection velocities of 0.025 m/s and 0.5 m/s the characteristic Reynolds numbers for pore body E are, respectively, about 12 and 180 in this case. Fig. 9 compares the computed velocity vector ﬁelds at the mid-section of pore body F, with one 1 mm, one 0.6 mm and two 0.2 mm connecting throats, for inlet velocities of 0.025 m/s and 0.5 m/s. The general features of the scaled velocity vector ﬁelds are roughly similar. The ﬂow enters the pore body from the lower-left 0.6 mm throat; and the bulk of the ﬂow leaves through the left 1 mm throat. The velocity magnitudes in the 0.2 mm throat are quite low. Fig. 9 also shows some distinct diﬀerences between the two ﬂow rates. In the lowvelocity case, the ﬂow turns smoothly into the 1 mm throat, with no recirculating ﬂow region. For the inlet velocity of 0.5 m/s, the high-momentum jet penetrates deeply into the pore body and reaches the opposite wall. In addition to the large recirculation region in the pore body, non-fully developed ﬂows are seen where the ﬂuid exits the pore body into each of the two smallest throats.

Fig. 8. Velocity vector ﬁeld at the mid-section of pore body (E) with two 1 mm, one 0.6 mm and one 0.2 mm connecting tubes for inlet velocities of 0.025 and 0.5 m/s.

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

1275

Fig. 9. Velocity vector ﬁeld at the mid-section of pore body (F) with two 0.2 mm one 0.6 and one 1 mm connecting tubes for inlet velocities of 0.025 and 0.5 m/s.

In this case, the corresponding Reynolds numbers for pore body F (based on the throat having the maximum velocity) are about 15 and 240 for the injection velocities of 0.025 m/s and 0.5 m/s, respectively. Comparing the characteristic Reynolds numbers in pore bodies E and F for diﬀerent injection velocities, it is seen that the Reynolds number increased roughly proportional to the injection rate; but the rate of increase is not quite linear. Examination of the results for diﬀerent pore bodies and diﬀerent ﬂow rates shows that the variations illustrated by Fig. 8 are typical, and in some case the ﬂow patterns show some diﬀerences similar to those illustrated in Fig. 9. For inlet velocities of 0.025 m/s and 0.5 m/s, Fig. 10 compares the velocity magnitude contours across the entire ﬂow-cell. The corresponding velocity vector ﬁelds in an enlarged segment of the ﬂow-cell are shown in Fig. 11. It is seen that the ﬂow in the ﬂow-cell is focused into rel-

atively few high-speed channels. Careful examination of these ﬁgures indicates that the high-speed channels are passages of lowest resistance formed by the largest available connecting throats. This observation of the ﬂow behavior is as expected, since the pressure drop is approximately proportional to the cube of the connecting tube diameter. That is, for the same ﬂow rate, a 0.2 mm throat oﬀers about nine times the resistance of a 0.6 mm throat and 25 times resistance of a 1 mm throat. Figs. 10 and 11 also show that the normalized ﬂow patterns for the 0.025 m/s and 0.5 m/s inlet velocities are quite similar. A careful examination of the details of the ﬂow shows that there are some apparently minor local diﬀerences between the high and low speed ﬂows such as those shown in Figs. 8 and 9, but the general features do not change. That is, the path of minimum resistance does change signiﬁcantly with the increase in the ﬂow rate, at least in the range of ﬂow rate studied.

Fig. 10. Comparison of the velocity magnitude contours for inlet velocities of 0.025 and 0.5 m/s.

1276

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

Fig. 11. Comparison of the velocity magnitude contours for inlet velocities of 0.025 and 0.5 m/s.

Fig. 12. Comparison of the histogram of normalized velocity for a range of ﬂow rates.

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

To provide a better understanding of the features of the ﬂow ﬁeld in the ﬂow-cell as a function of the ﬂow rate, histograms of velocity magnitudes in the cell were evaluated. Fig. 12 compares the histogram of the normalized velocity ﬁelds in the ﬂow-cell for injection velocities of 0.025, 0.075, 0.375, and 0.5 m/s. The histogram velocities are normalized by the respective injection velocities, and the frequencies are shown in logarithmic scale. Clearly, many of the normalized velocities are very close to zero, where the histograms have the highest peak. Fig. 12 also shows that there are secondary peaks at higher values of the normalized velocities. It is conjectured that these peaks correspond to the velocities in the channels of minimum resistance that formed in the cell. For all inlet velocities studied, the peaks of the local normalized velocities appear at about 0.2, 0.4, 0.7 and 1. These peaks are sharpened as the injection velocity increases. This observation indicates that the preferred paths do not change with the increase in the ﬂow rate, as also seen in Figs. 10 and 11.

8. Darcy and non-Darcy ﬂow regimes Figs. 8–11 illustrate in detail the ﬂow ﬁelds at the pore level for two diﬀerent injection rates. Although the ﬂow patterns were similar in many ways, diﬀerences (e.g. recirculation zones at the larger ﬂow velocity) were observed. Fig. 13 plots the variation of the pressure diﬀerence across the lattice-cell with the inlet velocity and the inlet Reynolds number. The inlet Reynolds number, ReI, was based on the inlet velocity and the diameter of the inlet tube. The corresponding Darcy Reynolds number, ReD, based on the Darcy velocity and the throat size is much smaller. Using a throat size of 0.6 mm, it follows that ReD 0:0016 ReI .

ð4Þ

Fig. 13 shows that the pressure diﬀerence increased linearly with the Darcy velocity (inlet velocity) for

1277

speeds less than 8.6 mm/s (0.125 m/s). For larger Darcy velocities, the pressure drop increase deviates from the linear variation and a correction that is proportional to the square of velocity is needed. At low injection velocities the entire cell was in a creeping ﬂow condition, so the pressure drop is linearly related to the velocity; but at higher velocities this condition no longer obtained. To ﬁt the simulation data for the pressure drop, the following equation is proposed: l q ð5Þ rP ¼ V þ V 2 k g Here V is the Darcy velocity, k is the permeability of the medium, and 1/g is the inertial factor (also called the non-Darcy ﬂow coeﬃcient). For k = 6.5 · 1010 m2 and g = 25 · 106 m2 the predictions of Eq. (5) are illustrated in Fig. 13 by a solid line. This ﬁgure shows that Eq. (5) predicts the simulation data with reasonable accuracy. For low Darcy velocities, the linear part of Eq. (5) given as l ð6Þ rP ¼ V k is also plotted in Fig. 13 by a dashed line. It appears that the linear equation can predict the pressure drop for Darcy velocities less than about 10 mm/s and the inlet Reynolds numbers less than about 100 (ReD < 016). Fig. 13 shows that for ﬂow through the lattice-cell with larger Reynolds number, the pressure drop-Darcy velocity relationship is non-linear. Comparing the results with those in Figs. 8 and 9, one ﬁnds that as the Reynolds number in pore bodies increased, deviation from Darcy ﬂow increased as the onset of recirculation zones appeared in more pore bodies. These results are similar to those reported by Pan et al. [25], which were obtained using the lattice Boltzmann method. It should also be emphasized here that due to the complexity of the geometry, Darcy behavior could be observed locally in certain regions of the lattice-cell with higher Reynolds numbers, but the overall cell behavior will deviate from Darcy condition. Figs. 8–11 show the formation of recirculation zones in some of the pore bodies and developing ﬂows at the exits of some pore bodies. Comparing Figs. 8–11 with Fig. 13, it is seen that the onset of formation of recirculation zones and non-fully developed ﬂows qualitatively coincides with the beginning of non-Darcy ﬂow. Thus, the simulation results relate the microscopic eﬀects to the macroscopic ﬂow behavior. 9. Conclusions

Fig. 13. Variation of pressure drop with Darcy velocity and inlet Reynolds number in the ﬂow-cell.

In this study, ﬂow through a lattice-cell was studied using a combination of computer simulations and laboratory experiments. Experimental data obtained by a micro-PIV (particle image velocimetry) technique served

1278

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279

as an experimental check on the accuracy of the simulations; the simulations allowed us to examine the eﬀects of diﬀerent conditions much more rapidly than would have been possible in the laboratory. On the basis of the results presented, the following conclusions are drawn: • The computer simulations of ﬂow through the latticecell captured many of the features of the ﬂow behavior. • While some discrepancies were noticed, the simulation results were in good agreement with the experimental data obtained by the micro-PIV technique. • The simulation results show that the ﬂow in the lattice-cell formed high-speed region in passages of lowest resistance. • Many features of the normalized ﬂow ﬁelds did not change appreciably with signiﬁcant increases in the ﬂow rate. However, the macroscopic change from Darcy to non-Darcy ﬂow could be detected at the microscopic level and was associated with the onset of recirculation zones in the pore bodies. • As the macroscopic ﬂow rate increased (Reynolds number in pore bodies increased), deviations from Darcy ﬂow increased as recirculation zones appeared in more and more bodies. Acknowledgments This work was supported by Fossil Energy, US Department of Energy. The support of the Ohio Coal Development Oﬃce (Columbus, Ohio) for the experiments is acknowledged by J.R. Kadambi and B.Z. Saylor. References [1] Bear J. Hydraulics of ground water. New York: McGraw-Hill Publ. Co.; 1979. [2] Carman PC. Fluid ﬂow through a granular bed. Trans Inst Chem Eng Lond 1937;15:150–6. [3] Chen S, Doolen GD. Lattice Boltzmann method for ﬂuid ﬂows. Annu Rev Fluid Mech 1998;30:329–64. [4] Chen J, Wilkinson D. Pore-scale viscous ﬁngering in porous media. Phys Rev Lett 1985;55:1892–5. [5] Collins RE. Flow of ﬂuids through porous materials. New York: Reinhold Publishing Corporation; 1961. [6] Darcy H. Les Foundaines Publiques de la Ville de Dijon. Paris: Dalmont; 1856. [7] Devasenathipathy S, Santiago JG, Wereley ST, Meinhart CD, Takehara K. Particle imagining techniques for microfabricated ﬂuidic systems. Exp Fluids 2003;34:504–14. [8] Dullien FAL. Porous media: ﬂuid transport and pore structure. 2nd ed. New York: Academic Press; 1992. [9] Ergun S. Fluid ﬂow through packed columns. Chem Eng Progr 1952;48:89–94. [10] Fenwick DH, Blunt MJ. Network modeling of three phase ﬂow in porous media. SPE J 1998;3:86–97.

[11] Ferer M, Gump J, Smith DH. Fractal nature of viscous ﬁngering in three-dimensional pore-level models. Phys Rev E 1996;53: 2502–8. [12] Ferer M, Bromhal GS, Smith DH. Pore-level modeling of immiscible drainage: validation in the inversion percolation and DLA limits. Physica A 2003;319:11–35. [13] Ferer M, Bromhal GS, Smith DH. Two-phase ﬂow in porous media: crossover from capillary ﬁngering to compact invasion. Phys Rev E 2005;71:026303. [14] Forchheimer P. Wasserbwegung durch Boden. Z Ver Deuts Ing 1901;45:1781–8. [15] van Genabeek O, Rothman DH. Macroscopic manifestations of microscopic ﬂows through porous media; phenomenology from simulation. Annu Rev Earth Planet Sci 1996;24:63–87. [16] Hughes RG, Blunt MJ. Pore-scale modeling of multiphase ﬂow in fractures and matrix/ fracture transfer. In: Society of Petroleum Engineering Conference, Houston, TX, October 3–6, SPE 56411, 1999. [17] Hui M-H, Blunt MJ. Pore-scale modeling of three-phase ﬂow and the eﬀects of wettability. In: Society of Petroleum Engineering Conference, Tulsa, OK, April 3–5, SPE 59309, 2000. [18] Koponen A, Kataja M, Timonen J. Permeability and eﬀective porosity of porous media. Phys Rev E 1997;56:3319–25. [19] Kozeny J. Ueber kapillare Leitung des Wessers im Boden. S-Ber Weiner Akad, Abt IIa 1927;136:271–306. [20] Lenormand R, Touboul E, Zarcone C. Numerical models and experiments on immiscible displacements in porous media. J Fluid Mech 1988;189:165–87. [21] Maier RS, Kroll DM, Kutsovsky YE, Davis HT, Bernard RS. Simulation of ﬂow through bead packs using the lattice Boltzmann method. Phys Fluids 1998;10(1):60–74. [22] Martys N, Chen H. Simulation of multicomponent ﬂuids in complex three-dimensional geometries by the Lattice Boltzmann method. Phys Rev E 1996;53:743–50. [23] Meinhart CD, Wereley ST, Santiago JG. PIV measurements of a microchannel ﬂow. Exp Fluids 1999;27:414–9. [24] Missbach A. Listy Cukrovar 1937;55:293. [25] Pan C, Hilpert M, Miller CT. Pore-scale modeling of saturated permeabilities in random sphere packings. Phys Rev E 2001;64: 066702. [26] Pan C, Prins JF, Miller CT. A high-performance lattice Boltzmann implementation to model ﬂow in porous media. Comput Phys Commun 2004;158:89–105. [27] Raﬀel M, Willert CE, Kompenhans J. Particle image velocimetry: A practical guide. Berlin, Heidelberg, New York: Springer; 1998. [28] Reeves PC. The development of pore-scale network models for the simulation of capillary pressure-saturation-interfacial area-relative permeability relationships in multi-ﬂuid porous media. PhD dissertation, Princeton University, Princeton, NJ, 1997. [29] Reeves PC, Celia MA. A Functional relationship between capillary pressure, saturation, and interfacial area as related by a pore-scale model. Water Resour Res 1997;32(8):2345–58. [30] Santiago JG, Wereley ST, Meinhart CD, Beebe DJ, Adrian RJ. Particle image velocimetry system for microﬂuidics. Exp Fluids 1998;25:316–9. [31] Scheidegger AE. Physics of ﬂow through porous media. Toronto: University of Toronto Press; 1974. [32] Singh M, Mohanty KK. Dynamic modeling of drainage through three-dimensional porous materials. Chem Eng Sci 2003;58(1): 1–18. [33] Smith DH, Ahmadi G, Ji BG, Ferer M. Experimental and numerical studies of gas-liquid displacement in ﬂow-cells, with application to carbon dioxide sequestration in brine ﬁelds. In: FEDSM2002-31296, Proceedings of ASME FEDSMÕ02, ASME Fluids Engineering Division Summer Meeting, Montreal, Quebec, Canada, July 14–18, 2002.

A.R. Mazaheri et al. / Advances in Water Resources 28 (2005) 1267–1279 [34] Succi S, Foti E, Higuera F. Simulation of three-dimensional ﬂow in porous media with the lattice Boltzmann method. Europhys Lett 1989;10:433–8. [35] Tsakiroglou CD. A methodology for the derivation of nonDarcian models for the ﬂow of generalized newtonian ﬂuids in porous media. J Non-Newtonian Fluid Mech 2002;105:79–110. [36] White AM. Pressure drops and loading velocities in packed towers. Trans Amer Inst Chem Eng 1935;31:390–408.

1279

[37] Zerai B, Saylor BZ, Kadambi JR, Oliver M, Mazaheri AR, Ahmadi G, et al. Flow characterization through a network cell using microparticle image velocimetry. Trans Porous Media 2005; 60:159–81. [38] Zhang D, Zhang R, Chen S, Soll WE. Pore scale study of ﬂow in porous media: Scale dependency, REV, and Statistical REV. Geophys Res Lett 2000;27(8):1195–8.