Computer simulation of flow through a lattice flow-cell model

Computer simulation of flow through a lattice flow-cell model

Advances in Water Resources 28 (2005) 1267–1279 www.elsevier.com/locate/advwatres Computer simulation of flow through a lattice flow-cell model A.R. Ma...

1MB Sizes 0 Downloads 31 Views

Advances in Water Resources 28 (2005) 1267–1279 www.elsevier.com/locate/advwatres

Computer simulation of flow through a lattice flow-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 flow through a network model of a porous medium, we report (1) solutions of the Navier–Stokes equation for the flow, (2) micro-particle imaging velocimetry (PIV) measurements of local flow 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-different (square cross-section) pore-throat sizes were randomly distributed throughout the network. An unstructured computational grid for flow through an identical network was developed and used to compute the local pressure gradients and flow vectors for several different (macroscopic) flow rates. Numerical solution results were compared with the experimental data, and good agreement was found. Cross-over from Darcy flow to inertial flow was observed in the computational results, and the permeability and inertia coefficients of the network were estimated. The development of inertial flow was seen as a ‘‘two-step’’ process: (1) recirculation zones appeared in more and more pore bodies as the flow rate was increased, and (2) the strengths of individual recirculation zones increased with flow 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 flows; Porous media

1. Introduction The long-term motivation of this study is to provide a better understanding of the fundamentals of single and multiphase flow 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 flows [1,5,8,31]. Studies of single-phase flow through permeable media usually begin with DarcyÕs phenomenological equation [6] of about 150 years ago. According to DarcyÕs law, the total flow through such a medium is proportional to both the pressure gradient along the macroscopic-average flow direction and the cross-sectional area of the medium. Details of the microscopic geometry or of the flow inside the medium, however, are not addressed. At higher flow velocities, the pressure gradient and flow 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 modified to include a second-order term in the velocity. Later, the idea of non-integer exponents for the relationship between pressure gradient and flow 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 define the Reynolds number in terms of pore diameter for porous media (whether natural or synthetic) that have no simple tubular pores. Others, exemplified by the work of Kozeny [19], Carman [2], and Ergun [9], have turned to theoretical descriptions of laminar and turbulent flow to elucidate the Forchheimer equation [31]. Lattice Boltzmann (LB) approaches currently offer perhaps the most popular and promising theoretical method for single-phase flow 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 flow 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 flow of viscous fluids. Although the LB method is a powerful technique, it is computationally demanding. Pan et al. [26] developed a new approach for modeling flow through porous media using the LB method with higher-performance, particularly, in floating 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 flow 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 flow of shear-thinning (non-Newtonian) fluids in porous media. In this work, the analytical models were then integrated into numerical simulations of single-phase flow through a two-dimensional, square-lattice, pore-network model. However, almost all of this literature treats problems of flow of two (or more) ‘‘phases,’’ that is, either two fluids that are miscible and thus have no capillary pressure, but have ‘‘not yet’’ mixed; or else two fluids that truly are two phases at chemical equilibrium with a well-

defined capillary pressure that is additive with the Poiseuille viscous pressure. It has been shown that such models can successfully produce the special two-phase flow cases of DLA (diffusion-limited aggregation) and invasion percolation (IP) [12], and the crossover from capillary fingering to compact invasion [13]. Such models have seldom (if ever) been used for flow of a singlefluid, but are capable of treating single-fluid flow as a ‘‘trivial’’ case. Detailed measurements for the velocity field have been of considerable interest for understanding of the complex features of fluid flows and for model verification. Particle image velocimetry (PIV) is now a wellestablished technique that provides the velocity vector field in a plane of flow. In this approach laser pulses illuminate a plane in the test section. The flow is seeded with tracer particles that reflect laser light, and the particle positions at each laser pulse are recorded and tracked for evaluating the instantaneous velocity field [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 flow through porous media is now made possible by a triad of recent developments: (1) network models, (2) computers and computational methods for finding non-analytic solutions of the Navier–Stokes equation, and (3) microPIV measurements. Briefly: (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 flow 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 different locations and flow rates. For the lattice Boltzmann approach, measurements of velocity vectors and/or solutions of the Navier– Stokes equation offer an experimental check on the computations and/or an opportunity to compare LB results with results from classical flow equations. For two-phase network models that assume Poiseuille flow, the computational and experimental approach of this paper suggest the possibility of ‘‘going back’’ to explore the limitations of the Poiseuille flow assumption for single-phase flow; such explorations could illuminate the limitations of the assumption for two-phase flow, and/or lead to extensions of this paperÕs approach to include two-phase flow through networks. This paper describes a complementary combination of computations and laboratory experiments on singlephase flows in a lattice network. Viscous flow 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 flow rates, and the details of the velocity and pressure fields were evaluated. The effective permeability of the cell was estimated and the cross-over from Darcy flow to inertial flow was examined. An identical experimental flow-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 flow-cell, and generally good agreement was found.

2. Flow-cell description The experimental 64 · 64 mm flow-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 flow entered the cell through a 2.5 mm tube into a small plenum that connected to the pore bodies. At the exit area, flow from the cell discharged into a small plenum and left the flow-cell through a 2.5 mm outlet tube. Fig. 1 shows a picture of the experimental flow-cell used in the study. Experimental and computational results are reported for pore bodies B, C, and D and the flow throats connected directly to each of these pore bodies. Each of these pore bodies was connected to four

Fig. 2. Schematic geometry of the flow-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 flow cells used by Ferer et al. [11–13] and Smith et al. [33], who performed flow simulations of a much different type.

3. Numerical simulations A three-dimensional computational grid for simulation of fluid flow through the flow-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 five 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 flow-cell. The Navier–Stokes and continuity equations governing the motion of an incompressible viscous fluid are given as   oui oui oP o2 ui q þ uj þl þ qgi ; ð1Þ ¼ oxi ot oxj oxj oxj

Fig. 1. Picture of the flow-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 flow-cell model.

Here, q is the density of the fluid, ui is the fluid velocity, P is the pressure, and l is the fluid viscosity. The FLUENTTM code (version 6.02) was used, and the flow 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 flow 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 difference 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 figure shows that the flow 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 flow-cell model.

40 mm/s and 15 mm/s. However, when the fluid leaves the pore body, through the two right-hand throats, the fluid 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 different throats are different, several different Reynolds numbers may be defined. We used the hydraulic diameter of the throat with maximum velocity and the corresponding average velocity in that throat to define 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 fluid. 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 flows in the throats are not necessarily fully developed with a well established parabolic profile. Fig. 5b shows the pressure contour plot across the pore body in the same plane as Fig. 5a. This figure shows that the pressure difference 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 flow 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 figure shows that the flow 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 flow 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 flow 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 difference between the pore body and the lower-right throat with negligible flow rate. This figure also shows a maximum pressure difference 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 figure shows a very low fluid velocity inside the pore body, while the liquid velocity is about 10 mm/s in all the connecting throats. This figure also shows that the flow 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 difference 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 fluid flow 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 flow 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 flow. Micro-particle image velocimetry (PIV) was used to measure the velocity fields 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 field of view (FOV), 2500 lm · 2500 lm and 5000 lm · 5000 lm, were used to measure the flow fields 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 profile in the 2.5 mm diameter test cell exit tube was laminar. The flow rate obtained from the measured velocity profile was within ±1.6% of the flow rate obtained by collecting and weighing the produced fluid. 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 flow 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 flow 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 flow direction).

Table 1 Laser, camera, fluid 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 different 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 figure shows that the general pattern of the measured velocity field is quite similar to the computed one shown in Fig. 5a. In particular, the direction and magnitude of the flow field 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 field 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

figure 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 inflow velocity of about 35 mm/s in the lower left inlet throat. The outflow 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 field inside the pore body. Nevertheless, the general features of the measured flow pattern are comparable with the computed one. In particular, as in the computer simulation results, the measured velocity field shows that the flow 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 flow-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 flow rates. The flow fields 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 fields. Here the results for inlet velocities of

0.025 m/s and 0.5 m/s are compared, so that the flow rate for the high speed case is 20 times the baseline flow rate. In Fig. 8, the velocity is scaled relative to its peak value at each flow rate, so that a comparative study can be performed. This figure shows that the general features of the normalized flow fields appear quite similar for the two flow rates. The flow 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 flow-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, different Reynolds number may be defined. Here we used Eq. (3) to define 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 fields 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 fields are roughly similar. The flow enters the pore body from the lower-left 0.6 mm throat; and the bulk of the flow 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 differences between the two flow rates. In the lowvelocity case, the flow turns smoothly into the 1 mm throat, with no recirculating flow 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 flows are seen where the fluid exits the pore body into each of the two smallest throats.

Fig. 8. Velocity vector field 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 field 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 different 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 different pore bodies and different flow rates shows that the variations illustrated by Fig. 8 are typical, and in some case the flow patterns show some differences 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 flow-cell. The corresponding velocity vector fields in an enlarged segment of the flow-cell are shown in Fig. 11. It is seen that the flow in the flow-cell is focused into rel-

atively few high-speed channels. Careful examination of these figures indicates that the high-speed channels are passages of lowest resistance formed by the largest available connecting throats. This observation of the flow behavior is as expected, since the pressure drop is approximately proportional to the cube of the connecting tube diameter. That is, for the same flow rate, a 0.2 mm throat offers 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 flow 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 flow shows that there are some apparently minor local differences between the high and low speed flows 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 significantly with the increase in the flow rate, at least in the range of flow 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 flow rates.

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

To provide a better understanding of the features of the flow field in the flow-cell as a function of the flow rate, histograms of velocity magnitudes in the cell were evaluated. Fig. 12 compares the histogram of the normalized velocity fields in the flow-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 flow rate, as also seen in Figs. 10 and 11.

8. Darcy and non-Darcy flow regimes Figs. 8–11 illustrate in detail the flow fields at the pore level for two different injection rates. Although the flow patterns were similar in many ways, differences (e.g. recirculation zones at the larger flow velocity) were observed. Fig. 13 plots the variation of the pressure difference 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 difference 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 flow condition, so the pressure drop is linearly related to the velocity; but at higher velocities this condition no longer obtained. To fit 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 flow coefficient). 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 figure 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 flow 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 finds that as the Reynolds number in pore bodies increased, deviation from Darcy flow 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 flows 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 flows qualitatively coincides with the beginning of non-Darcy flow. Thus, the simulation results relate the microscopic effects to the macroscopic flow behavior. 9. Conclusions

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

In this study, flow 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 effects of different 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 flow through the latticecell captured many of the features of the flow 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 flow in the lattice-cell formed high-speed region in passages of lowest resistance. • Many features of the normalized flow fields did not change appreciably with significant increases in the flow rate. However, the macroscopic change from Darcy to non-Darcy flow could be detected at the microscopic level and was associated with the onset of recirculation zones in the pore bodies. • As the macroscopic flow rate increased (Reynolds number in pore bodies increased), deviations from Darcy flow 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 Office (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 flow through a granular bed. Trans Inst Chem Eng Lond 1937;15:150–6. [3] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 1998;30:329–64. [4] Chen J, Wilkinson D. Pore-scale viscous fingering in porous media. Phys Rev Lett 1985;55:1892–5. [5] Collins RE. Flow of fluids 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 fluidic systems. Exp Fluids 2003;34:504–14. [8] Dullien FAL. Porous media: fluid transport and pore structure. 2nd ed. New York: Academic Press; 1992. [9] Ergun S. Fluid flow through packed columns. Chem Eng Progr 1952;48:89–94. [10] Fenwick DH, Blunt MJ. Network modeling of three phase flow in porous media. SPE J 1998;3:86–97.

[11] Ferer M, Gump J, Smith DH. Fractal nature of viscous fingering 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 flow in porous media: crossover from capillary fingering 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 flows 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 flow 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 flow and the effects 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 effective 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 flow through bead packs using the lattice Boltzmann method. Phys Fluids 1998;10(1):60–74. [22] Martys N, Chen H. Simulation of multicomponent fluids 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 flow. 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 flow in porous media. Comput Phys Commun 2004;158:89–105. [27] Raffel 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-fluid 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 microfluidics. Exp Fluids 1998;25:316–9. [31] Scheidegger AE. Physics of flow 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 flow-cells, with application to carbon dioxide sequestration in brine fields. 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 flow 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 flow of generalized newtonian fluids 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 flow in porous media: Scale dependency, REV, and Statistical REV. Geophys Res Lett 2000;27(8):1195–8.