PHYSICAL REVIEW E 88, 063014 (2013)
Hyperbolic regions in flows through three-dimensional pore structures Jeffrey D. Hyman* University of Arizona, Program in Applied Mathematics, Tucson, Arizona 85721-0089, USA
C. Larrabee Winter† University of Arizona, Department of Hydrology and Water Resources, Program in Applied Mathematics, Tucson, Arizona 85721-0011, USA (Received 6 June 2013; published 19 December 2013) Finite time Lyapunov exponents are used to determine expanding, contracting, and hyperbolic regions in computational simulations of laminar steady-state fluid flows within realistic three dimensional pore structures embedded within an impermeable matrix. These regions correspond approximately to pores where flow converges (contraction) or diverges (expansion), and to throats between pores where the flow mixes (hyperbolic). The regions are sparse and disjoint from one another, occupying only a small percentage of the pore space. Nonetheless, nearly every percolating fluid particle trajectory passes through several hyperbolic regions indicating that the effects of in-pore mixing are distributed throughout an entire pore structure. Furthermore, the observed range of fluid dynamics evidences two scales of heterogeneity within each of these flow fields. There is a larger scale that affects dispersion of fluid particle trajectories across the connected network of pores and a relatively small scale of nonuniform distributions of velocities within an individual pore. DOI: 10.1103/PhysRevE.88.063014
PACS number(s): 47.56.+r, 45.70.−n, 47.51.+a, 47.52.+j
I. INTRODUCTION
The study of fluid flow through porous media can be roughly partitioned into two classes based on the associated length scales. A larger (Darcy) scale (cm-km) represents flow as occurring through a continuum characterized by a parameter, hydraulic conductivity, that is effectively the inverse of the medium’s resistance to flow. There are also flows through porous media represented as explicit structures of interconnected pores where there are two even smaller length scales. A mesoscopic scale (mm) corresponds to the effects of the geometry and topology of the pore structure on the trajectories of fluid particles as they transit the pore space. A smaller, microscopic, scale (microns) is needed to capture the effects of heterogeneities of flow fields occurring within individual pores. The effects of both scales have been recognized for some time [1,2], but their relative effects have been difficult to observe or analyze. In contrast to the Darcy scale where theoretical, computational, and experimental techniques are plentiful, flow through interconnected pore structures is difficult to solve for analytically due to the complexity of the pore geometry and topology or to observe through physical experiments due to material constraints. Until recently the application of computational fluid dynamics in this area was limited to relatively simple representations of the pore structure, e.g., regular lattices or pore volumes [3–8], despite the relevance of flow through realistic pore structures to many areas of science and engineering including geophysical processes, filter design, petroleum extraction, nuclear waste disposal, and carbon sequestration. Recent computational advances are allowing researchers to link various flow properties to the physical attributes of more complex and realistic pore structures [9–13], but the influence of the pore geometry and network topology
* †
Corresponding author: [email protected] [email protected]
1539-3755/2013/88(6)/063014(6)
on fluid dynamics and the mixing of solutes therein is not well understood. Models of fluid flow through explicit networks of pores frequently assume that mixing occurs uniformly at every scale and throughout the pore space, even though observations have shown that this conception is incorrect when the effects of diffusion are negligible, e.g., high P´eclet numbers or a pre-asymptotic regime [1,14,15]. Computational experiments reported here indicate that mixing is induced in such cases by mechanical dispersion, which is a function of both the microscopic scale of in-pore fluid dynamics and the mesoscopic scale of entire pore structure. In particular, variations in the microscopic local pore geometry and mesoscopic network topology interact to induce nonuniform fluid flow fields characterized by regions of confluence and divergence. Regions of contraction and expansion in the flow field correspond to stable and unstable manifolds in the velocity field, respectively, and their intersections are hyperbolic regions where mixing occurs [16,17]. Although the boundaries of the solid matrix are generally too complicated to extract the precise locations of these manifolds, finite time Lyapunov exponents (FTLE) can be used to identify corresponding regions of expansion and contraction [18]. In this paper, FTLE fields of steady-state three dimensional laminar flows through four independent simulations of realistic porous media with varying permeabilities are computed and their relationship with the flow field and pore structure analyzed. The forward and backward FTLE fields are used to decompose each velocity field into regions of expansion, contraction, both, or neither thereby identifying confluence, divergence, and mixing regions in the flow. Regions of high values in the forward and backward FTLE fields are sparsely distributed throughout the pore spaces, and their locations are found to be correlated with larger velocities in the primary direction of flow and the absolute value of horizontal divergence which indicates lateral movement. Even sparser are hyperbolic regions identified by the intersections of the expanding and contracting regions. Hyperbolic regions
063014-1
©2013 American Physical Society
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013)
are found in pores where junctions in the network induce both confluence and divergence in the same place. Less than ten percent of the void volumes are occupied by a hyperbolic region suggesting that advective mixing in flows through complicated pore structures is isolated. Nonetheless, greater than 99% of all percolating fluid particle trajectories intersect at least one hyperbolic region indicating that the effects of localized mixing are distributed throughout the entire pore space.
TABLE I. Characteristics of sample pore spaces: sample number (Si ), porosity (n), mean hydraulic radius (R; m), tortuosity (τ ), permeability (k; m2 ).
n R τ k
S1
S2
S3
S4
0.29 6.64 × 10−5 1.38 4.12 × 10−10
0.33 7.10 × 10−5 1.36 6.12 × 10−10
0.38 7.77 × 10−5 1.33 8.73 × 10−10
0.45 8.75 × 10−5 1.27 1.41 × 10−9
II. PORE SPACE GENERATION AND SIMULATION OF FLOW
III. FINITE TIME LYAPUNOV EXPONENTS
The methods of Smolarkiewicz and Winter [19] and Hyman et al. [11] are used to simulate steady-state laminar fluid flow within stochastically generated isotropic and statistically stationary three dimensional porous media. To create a pore space realization, first a correlated random topography is produced on a grid by convolving a field of uniform random variables with a Gaussian kernel. Then a level threshold is applied to the resulting topography to determine which locations are in the pore space and which are in the solid matrix. Each three dimensional virtual porous medium has dimensions Lx = Ly = Lz = 1.27 × 10−2 m, with volume V = 2.05 × 10−6 m3 . The domain is discretized using a uniform Cartesian grid with a grid spacing of δx = 10−4 m. Each medium is periodic in the vertical direction with no flow boundary conditions applied along lateral boundaries. The Navier-Stokes equations, ∇ · v = 0, ∂v + v · ∇v = −∇π + g + μv − αv, ∂t
(1)
are numerically integrated using the EULAG system [20] to determine the fluid velocity field, v(x), within the void space of a porous medium. The primes in (1) refer to perturbations with respect to static ambient atmospheric conditions characterized by a constant density ρ0 and pressure p0 = p0 (z), so π = (p − p0 )/ρ and g = (0,0, − gρ /ρ) where ρ = const ρ0 denotes the density of fluid and g is the gravitational acceleration, g = 9.81 m/s2 . The kinematic viscosity of water μ is 10−6 ms−2 . No slip internal boundary conditions along pore walls are enforced using an immersed boundary method [21] by inserting a fictitious body force into the equations of motion; α in (1). The resulting dynamics are such that velocity is negligible and pressure irrelevant within the solid matrix where the body forces are high. The flow simulations are run for 5 × 10−2 sec with steady state conditions reached in 2–3 × 10−2 sec. Characteristics of the four sample pore spaces, Si , are computed using the methods of Hyman et al. [22] and are reported in Table I. The porosities fall at the lower end of what is considered a normal range, 0.2 n 0.7, and the mean hydraulic radius, tortuosities, and permeabilities correspond to well-sorted gravel or sand [22]. The measurements of tortuosity reported in Table I indicate that fluid particle trajectories deviate from a straight line path, for which τ = 1.
The scalar field of finite time Lyapunov exponents (FTLE) is computed everywhere in the pore space according to the techniques of Shadden et al. [18]. The FTLE, σxT0 , at a given point, x0 , measures the degree of stretching experienced by an infinitesimal fluid element advected with the flow for a finite time T ; larger values indicate more stretching. Specifically, σxT0 = where
1 ln λmax (S(x0 ,T )), 2|T |
S(x0 ,T ) =
dφ T (x(t)) dx(t)
∗
dφ T (x(t)) dx(t)
(2) (3)
is the (right) Cauchy-Green strain tensor, λmax is the maximum eigenvalue of S(x0 ,T ), and A∗ denotes the adjoint of the matrix A. The absolute value is used in (2) to allow for both positive and negative values of T . The strain tensor associated with x0 and T is constructed using the gradient of the flow map φ, which is determined by numerically integrating the trajectory equation, x˙ (x0 ,t) = v(x(x0 ,t)),
(4)
with a fourth order Runge-Kutta scheme coupled with trilinear interpolation to compute the velocities away from computational nodes. The sample pore volume S2 is used to determine a suitable value of T so that σxT0 values have stabilized while still capturing local dynamics, T = 2 × 10−2 sec. The forward flow map, T > 0, is used to identify expansion, while the backward flow map, T < 0, is used to identify contraction. The left side of Fig. 1 is a contour plot of the forward (expanding) FTLE field and the right side is a contour plot of the backward (contracting) FTLE field from one-fourth of a horizontal cross section of sample pore volume S3 (cf. Table I). Contour plots from the other samples are similar. High values in the forward FTLE field form disjoint clusters, and similar characteristics are observed in the backward FTLE field, but those clusters are sparser than clusters observed in the forward FTLE field. This fragmentation indicates that significant expansion and contraction occur in isolated regions within the pore space, and demonstrates the nonuniformity of the flow within pores and across the pore space. To identify locations where significant expansion or contraction in the flow field occur, the forward and backward FTLE fields are thresholded according to σxT0 σ¯ xT + C σˆ xT − σ¯ xT , (5)
063014-2
HYPERBOLIC REGIONS IN FLOWS THROUGH THREE- . . .
PHYSICAL REVIEW E 88, 063014 (2013)
1mm 0
2
4
6
8
10
12
14
16
−3
18 x 10
FIG. 1. (Color online) Contour plots of the forward (left) and backward (right) FTLE fields in one-fourth of a horizontal cross section from a porous medium with porosity 0.38 show that regions of high FTLE values are fragmented. Solid matrix shown in black.
which partitions each field into high and low values [23]. Here, σ¯ xT and σˆ xT are the mean and maximum FTLE values, respectively, and the constant C determines the minimum FTLE value that constitutes significant expansion or contraction. A region of expansion is a connected cluster composed of at least 50 computational grid nodes in the forward FTLE field whose values are greater than the threshold level. Regions of contraction are defined with the same criteria applied to the backward FTLE field. The sample pore volume S2 was used to determine that a value of C = 0.10 extracts the expanding and contracting regions from the respective FTLE fields well. Restricting analysis to clusters composed of at least 50 nodes removes negligible volumes. Table II reports the percentage of the void volume occupied by expanding and contracting regions, the number of regions in each sample, and the correlation coefficient of the expanding and contracting regions with vertical velocity, the primary direction of flow, and the absolute value of horizontal divergence. Expanding regions occupy a larger percentage of the void space than contracting regions indicating that fluid particle trajectories are more prone to disperse than coalesce, and that percentage rises with permeability, porosity, and hydraulic radius. Combined with the decline in the number of regions, this increase indicates that the volume of the expanding regions increases with permeability; images confirm this growth. Across all of the samples, the expanding regions are more
correlated with the vertical velocity and horizontal divergence fields than the contracting regions. In all four samples, the percentage of the void space occupied by the contracting regions is small, less than 20%. The sample with the lowest permeability, S1 , has nearly double the percentage of its void space occupied by these regions compared to the other three samples, indicating that the flow coalesces into select channels through this pore network. This result is consistent with the findings of Hyman et al. [11] who observed that a high percentage of the flux through horizontal cross sections occurred in a low percentage of the void space in low porosity samples. Although the percentage of the void space occupied by the contracting regions remains roughly constant in the higher permeability samples, the number of regions increases. Hence, the regions are smaller in higher permeability samples; images confirm this decrease in volume. The regions of contraction are more correlated with both the vertical velocity and horizontal divergence in lower permeability porous media than in samples with higher permeabilities. The low correlations between the vertical velocity and the absolute value of the horizontal divergence with the expanding and contracting regions suggest that neither of these two quantities of the velocity field alone can fully account for the measured confluence and/or divergence throughout the pore structure. However, because the flow is incompressible, the flow map is volume preserving and the flow field must account for all expansion and contraction. One explanation for these low correlations is that the underlying cause of expansion or contraction depends upon different combinations of the horizontal divergence and vertical velocity and varies spatially. The disparity between locations results from the variable resistance offered by the mesoscopic geometry and topology of the pore structure. The influence of variations in the pore structure on the FTLE fields is demonstrated in Fig. 2 which shows a closeup of the lower right corner of the thresholded forward FTLE field extracted from the left cross section in Fig. 1 at four sequential vertical levels, descending clockwise from the upper left corner. As the height decreases, the pore splits into three and a region of expansion is observed at this junction. Similar behavior is observed in the backward FTLE field where regions of contraction coincide with a convergence in the flow that occurs at pore junctions. IV. HYPERBOLIC REGIONS
TABLE II. Statistics of expanding and contracting regions: sample number (Si ), percent of void volume occupied (PV ), number of regions (NR ), correlation coefficient with vertical velocity (rw ), correlation coefficient with the absolute value of horizontal divergence (rH ). Expansion
Contraction
Si
PV
NR
rw
rH
PV
NR
rw
rH
S1 S2 S3 S4
24.82 27.06 27.76 28.44
176 155 134 137
−0.37 −0.34 −0.31 −0.28
0.36 0.33 0.28 0.23
18.29 8.79 12.24 9.59
251 307 381 389
−0.35 −0.27 −0.28 −0.24
0.33 0.22 0.20 0.12
Hyperbolic regions occur where stable and unstable manifolds in the flow field intersect and are identified by locations that are in both expanding and contracting regions. To provide intuition about the pore geometry responsible for these hyperbolic regions, the left side of Fig. 3 shows a subvolume (about one-thousandth the entire volume) of the pore volume S4 where fluid particle trajectories from two separate pores (red and blue indicate different initial pores) converge into a hyperbolic region (yellow spheres) and then disperse back into the pore network; the primary direction of flow is from top to bottom. The trajectories are separated before they enter the hyperbolic region, intermingle therein, and are mixed together when they disperse into two disjoint pores below.
063014-3
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013) TABLE III. Statistics of hyperbolic regions; sample number (Si ), percent of void volume occupied (PV ), percent of percolating particles that intersect at least one hyperbolic regions (PT ), averages given a particle intersects a hyperbolic region: number of unique intersections (NT ), percent of particle travel time in a hyperbolic region (TT ).
−3
x 10 18 16 14
Unrestricted
Restricted
12
Si
PV
PT
NT
TT
PV
PT
NT
TT
10
S1 S2 S3 S4
9.46 4.81 5.87 4.50
100.00 99.11 100.00 99.07
12.96 10.63 12.06 8.49
17.50 7.07 7.57 5.90
4.43 1.05 1.12 0.61
100.00 65.18 69.49 51.01
3.43 1.88 1.95 1.31
10.64 2.54 2.65 1.98
8 6
this phenomena of in-pore mixing requires that the flow under consideration be three dimensional. In a similar study involving simpler geometry and topology, Cardenas [9] observed vortices near pore throats when considering three dimensional Navier-Stokes flow through spherical and ellipsoidal cubicpacked grains at a range of Reynolds numbers. Table III reports statistics of hyperbolic regions and their intersections with fluid particle trajectories that percolate through the entire domain, determined using the techniques of Hyman et al. [11]. Unrestricted values are generated using all hyperbolic regions, and restricted values are the hyperbolic regions which are composed of at least 50 nodes. The percentages of the void volume occupied by the hyperbolic regions are less than the corresponding values of either the expanding or contracting regions (cf. PV Table II), decrease with increasing permeability, and demonstrate the sparsity of hyperbolic regions, especially those with larger volumes as evidenced by the even smaller percentages reported for the restricted regions. Despite the isolation of the hyperbolic regions that results from the fragmentation of the expanding and contracting regions, they are connected via fluid particle trajectories. Hence, the product of this in-pore mixing is distributed throughout the rest of the pore space via streamlines in the flow. Even though the hyperbolic regions occupy less than
4 2 0
1 mm FIG. 2. (Color online) Closeup of the lower right corner of the thresholded forward FTLE field extracted from the left cross section in Fig. 1 at four sequential vertical levels, descending clockwise from the upper left corner. The thresholded forward FTLE field shows that a region of expansion occurs where a pore splits and the flow separates.
The subfigure on the right is a closeup of the hyperbolic regions on the left with selected trajectories colored yellow. Starting from the black cross in the center of the image, placed at the center of the hyperbolic region, the paths of the yellow trajectories are integrated backwards in time to show where they originate and forward in time to show where they go. The yellow trajectories originate in separate pores above and disperse into two disjoint pores below, thereby demonstrating the pore-scale mixing which occurs in hyperbolic regions. Furthermore, because the flow is laminar and at steady state
1 mm FIG. 3. (Color online) Left: Subvolume of the sample S4 (about one-thousandth the entire volume) with fluid particle trajectories from two separate pores (red and blue indicate different initial pores) converging into a hyperbolic region (yellow spheres) and diverging back into the pore network. Right: Closeup view of the of hyperbolic region on the left. Starting from the black cross in the center of the image, placed at the center of the hyperbolic region, the paths of the yellow trajectories are integrated backwards in time to show where they originate and forward in time to show where they exit the subvolume. 063014-4
HYPERBOLIC REGIONS IN FLOWS THROUGH THREE- . . .
PHYSICAL REVIEW E 88, 063014 (2013)
10% of the void volume, greater than 99% of percolating fluid particles pass through a control volume whose corresponding grid node is in a hyperbolic region. Those that do not remain close to the impermeable lateral boundaries of the domain throughout their trajectories. The lower values reported for the restricted regions indicate that intersections with larger hyperbolic regions are less frequent. Nonetheless, at least half of the particles intersect between one and two of these larger hyperbolic regions, and in the lowest permeability sample all percolating particles still pass through at least one of the larger hyperbolic regions. Moreover, particles that intersect a hyperbolic region pass through ten unique regions on average and spend a nontrivial percent of their travel time therein. Additional statistics not given show that particles intersect a hyperbolic region shortly after insertion. The probability that a particle’s first intersection with a hyperbolic region occurs at a given height decays exponentially with increasing vertical distance from the particle’s insertion point. The first-passage time distributions reported by Hyman et al. [11] illustrate that flows through virtual pore structures sampled from the same ensemble as those used here demonstrate longitudinal dispersion whose magnitude decreases with porosity. The measurements of tortuosity reported in Table I indicate that fluid particle trajectories deviate from a straight line path and indicate that transverse dispersion is also occurring. Computation of particles’ lateral displacement from their initial position shows that on average, particles disperse horizontally at least one-tenth of the domain width as they pass from the top of the domain to the bottom in all of the samples considered. In conjunction with the high percentage of particle intersections with hyperbolic regions (cf. PT in Table III) this suggests that the effects of localized mixing are dispersed laterally as well as in the primary direction of flow through an entire pore structure.
flows through realistic three dimensional pore structures are sparse and disjoint. They coincide with pore throats where both confluence and divergence in the flow occur, and occupy a low percentage of the void volume. Nonetheless, nearly every percolating fluid particle trajectory passes through several hyperbolic regions and together their trajectories form pathways between these mixing regions. The percolating trajectories and hyperbolic regions provide a skeleton of mixing in the flow field and indicate that the effects of localized mixing are distributed across trajectories and throughout an entire pore space. Figure 2 shows a region of expansion occurring where a pore splits, and highlights how the resistance offered by the intricately connected pore structure induces a nonuniform flow therein. Through Lagrangian analysis it was determined that this nonuniformity occurs on two length scales in the flow field. There is a larger scale of nonuniformity across the entire flow field that corresponds to what is classically refereed to as mechanical dispersion. This is evidenced by the contour plots of the forward and backward FTLE fields (Fig. 1), tortuosities of percolating particles (Table I), and their lateral displacement (Sec. IV). There is also a relatively small scale that occurs in individual pores, evidenced by the mixing observed in localized hyperbolic regions (Fig. 3) that has not been extensively investigated previously due to limitations in the ability to observe flows within pore spaces on the one hand, and computational limitations on the other. Advective tracer particles that percolate through the entire domain link the two lengths scales together (Table III) indicating that mixing and transport results from the interaction of these two scales of dispersion.
V. SUMMARY AND REMARKS
ACKNOWLEDGMENTS
Through analysis of finite time Lyapunov exponent fields it was determined that hyperbolic regions in steady state laminar
We would like to thank R. Indik, K. Lin, and M. Tabor for discussions and insightful comments.
[1] J. Bear, Dynamics of Fluids in Porous Media (Dover, New York, 1988). [2] M. Sahimi, Flow and Transport in Porous Media and Fractured Rock (Wiley, Weinheim, 2012). [3] H. Hasimoto, J. Fluid Mech. 5, 317 (1959). [4] B. E. Launder and T. H. Massey, J. Heat Transf. 100, 565 (1978). [5] A. S. Sangani and A. Acrivos, Int. J. Multiphase Flow 8, 343 (1982). [6] A. S. Sangani and A. Acrivos, Int. J. Multiphase Flow 8, 193 (1982). [7] A. Zick and G. M. Homsy, J. Fluid Mech. 115, 13 (1982). [8] A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 54, 406 (1996). [9] M. B. Cardenas, Geophys. Res. Lett. 35, L18402 (2008).
[10] B. Bijeljic, P. Mostaghimi, and M. J. Blunt, Phys. Rev. Lett. 107, 204502 (2011). [11] J. D. Hyman, P. K. Smolarkiewicz, and C. L. Winter, Phys. Rev. E 86, 056701 (2012). [12] B. Bijeljic, A. Raeini, P. Mostaghimi, and M. J. Blunt, Phys. Rev. E 87, 013011 (2013). [13] M. Matyka, Z. Koza, J. Gołembiewski, M. Kostur, and M. Januszewski, Phys. Rev. E 88, 023018 (2013). [14] S. Neuman, Adv. Water Resour. 28, 149 (2005). [15] A. M. Tartakovsky and S. P. Neuman, Geophys. Res. Lett. 35, L21401 (2008). [16] J. M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport (Cambridge University Press, Cambridge, England, 1989), Vol. 3. [17] G. A. Voth, G. Haller, and J. P. Gollub, Phys. Rev. Lett. 88, 254501 (2002).
063014-5
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013)
[18] S. C. Shadden, F. Lekien, and J. E. Marsden, Physica D 212, 271 (2005). [19] P. K. Smolarkiewicz and C. L. Winter, J. Comput. Phys. 229, 3121 (2010). [20] J. M. Prusa, P. K. Smolarkiewicz, and A. A. Wyszogrodzki, Comput. Fluids 37, 1193 (2008).
[21] R. Mittal and G. Iaccarino, Annu. Rev. Fluid Mech. 37, 239 (2005). [22] J. D. Hyman, P. K. Smolarkiewicz, and C. L. Winter, Water Resour. Res. 49, 2080 (2013). [23] Y. Zhu and R. Newell, Mon. Weather Rev. 126, 725 (1998).
063014-6
Hyperbolic regions in flows through three-dimensional pore structures Jeffrey D. Hyman* University of Arizona, Program in Applied Mathematics, Tucson, Arizona 85721-0089, USA
C. Larrabee Winter† University of Arizona, Department of Hydrology and Water Resources, Program in Applied Mathematics, Tucson, Arizona 85721-0011, USA (Received 6 June 2013; published 19 December 2013) Finite time Lyapunov exponents are used to determine expanding, contracting, and hyperbolic regions in computational simulations of laminar steady-state fluid flows within realistic three dimensional pore structures embedded within an impermeable matrix. These regions correspond approximately to pores where flow converges (contraction) or diverges (expansion), and to throats between pores where the flow mixes (hyperbolic). The regions are sparse and disjoint from one another, occupying only a small percentage of the pore space. Nonetheless, nearly every percolating fluid particle trajectory passes through several hyperbolic regions indicating that the effects of in-pore mixing are distributed throughout an entire pore structure. Furthermore, the observed range of fluid dynamics evidences two scales of heterogeneity within each of these flow fields. There is a larger scale that affects dispersion of fluid particle trajectories across the connected network of pores and a relatively small scale of nonuniform distributions of velocities within an individual pore. DOI: 10.1103/PhysRevE.88.063014
PACS number(s): 47.56.+r, 45.70.−n, 47.51.+a, 47.52.+j
I. INTRODUCTION
The study of fluid flow through porous media can be roughly partitioned into two classes based on the associated length scales. A larger (Darcy) scale (cm-km) represents flow as occurring through a continuum characterized by a parameter, hydraulic conductivity, that is effectively the inverse of the medium’s resistance to flow. There are also flows through porous media represented as explicit structures of interconnected pores where there are two even smaller length scales. A mesoscopic scale (mm) corresponds to the effects of the geometry and topology of the pore structure on the trajectories of fluid particles as they transit the pore space. A smaller, microscopic, scale (microns) is needed to capture the effects of heterogeneities of flow fields occurring within individual pores. The effects of both scales have been recognized for some time [1,2], but their relative effects have been difficult to observe or analyze. In contrast to the Darcy scale where theoretical, computational, and experimental techniques are plentiful, flow through interconnected pore structures is difficult to solve for analytically due to the complexity of the pore geometry and topology or to observe through physical experiments due to material constraints. Until recently the application of computational fluid dynamics in this area was limited to relatively simple representations of the pore structure, e.g., regular lattices or pore volumes [3–8], despite the relevance of flow through realistic pore structures to many areas of science and engineering including geophysical processes, filter design, petroleum extraction, nuclear waste disposal, and carbon sequestration. Recent computational advances are allowing researchers to link various flow properties to the physical attributes of more complex and realistic pore structures [9–13], but the influence of the pore geometry and network topology
* †
Corresponding author: [email protected] [email protected]
1539-3755/2013/88(6)/063014(6)
on fluid dynamics and the mixing of solutes therein is not well understood. Models of fluid flow through explicit networks of pores frequently assume that mixing occurs uniformly at every scale and throughout the pore space, even though observations have shown that this conception is incorrect when the effects of diffusion are negligible, e.g., high P´eclet numbers or a pre-asymptotic regime [1,14,15]. Computational experiments reported here indicate that mixing is induced in such cases by mechanical dispersion, which is a function of both the microscopic scale of in-pore fluid dynamics and the mesoscopic scale of entire pore structure. In particular, variations in the microscopic local pore geometry and mesoscopic network topology interact to induce nonuniform fluid flow fields characterized by regions of confluence and divergence. Regions of contraction and expansion in the flow field correspond to stable and unstable manifolds in the velocity field, respectively, and their intersections are hyperbolic regions where mixing occurs [16,17]. Although the boundaries of the solid matrix are generally too complicated to extract the precise locations of these manifolds, finite time Lyapunov exponents (FTLE) can be used to identify corresponding regions of expansion and contraction [18]. In this paper, FTLE fields of steady-state three dimensional laminar flows through four independent simulations of realistic porous media with varying permeabilities are computed and their relationship with the flow field and pore structure analyzed. The forward and backward FTLE fields are used to decompose each velocity field into regions of expansion, contraction, both, or neither thereby identifying confluence, divergence, and mixing regions in the flow. Regions of high values in the forward and backward FTLE fields are sparsely distributed throughout the pore spaces, and their locations are found to be correlated with larger velocities in the primary direction of flow and the absolute value of horizontal divergence which indicates lateral movement. Even sparser are hyperbolic regions identified by the intersections of the expanding and contracting regions. Hyperbolic regions
063014-1
©2013 American Physical Society
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013)
are found in pores where junctions in the network induce both confluence and divergence in the same place. Less than ten percent of the void volumes are occupied by a hyperbolic region suggesting that advective mixing in flows through complicated pore structures is isolated. Nonetheless, greater than 99% of all percolating fluid particle trajectories intersect at least one hyperbolic region indicating that the effects of localized mixing are distributed throughout the entire pore space.
TABLE I. Characteristics of sample pore spaces: sample number (Si ), porosity (n), mean hydraulic radius (R; m), tortuosity (τ ), permeability (k; m2 ).
n R τ k
S1
S2
S3
S4
0.29 6.64 × 10−5 1.38 4.12 × 10−10
0.33 7.10 × 10−5 1.36 6.12 × 10−10
0.38 7.77 × 10−5 1.33 8.73 × 10−10
0.45 8.75 × 10−5 1.27 1.41 × 10−9
II. PORE SPACE GENERATION AND SIMULATION OF FLOW
III. FINITE TIME LYAPUNOV EXPONENTS
The methods of Smolarkiewicz and Winter [19] and Hyman et al. [11] are used to simulate steady-state laminar fluid flow within stochastically generated isotropic and statistically stationary three dimensional porous media. To create a pore space realization, first a correlated random topography is produced on a grid by convolving a field of uniform random variables with a Gaussian kernel. Then a level threshold is applied to the resulting topography to determine which locations are in the pore space and which are in the solid matrix. Each three dimensional virtual porous medium has dimensions Lx = Ly = Lz = 1.27 × 10−2 m, with volume V = 2.05 × 10−6 m3 . The domain is discretized using a uniform Cartesian grid with a grid spacing of δx = 10−4 m. Each medium is periodic in the vertical direction with no flow boundary conditions applied along lateral boundaries. The Navier-Stokes equations, ∇ · v = 0, ∂v + v · ∇v = −∇π + g + μv − αv, ∂t
(1)
are numerically integrated using the EULAG system [20] to determine the fluid velocity field, v(x), within the void space of a porous medium. The primes in (1) refer to perturbations with respect to static ambient atmospheric conditions characterized by a constant density ρ0 and pressure p0 = p0 (z), so π = (p − p0 )/ρ and g = (0,0, − gρ /ρ) where ρ = const ρ0 denotes the density of fluid and g is the gravitational acceleration, g = 9.81 m/s2 . The kinematic viscosity of water μ is 10−6 ms−2 . No slip internal boundary conditions along pore walls are enforced using an immersed boundary method [21] by inserting a fictitious body force into the equations of motion; α in (1). The resulting dynamics are such that velocity is negligible and pressure irrelevant within the solid matrix where the body forces are high. The flow simulations are run for 5 × 10−2 sec with steady state conditions reached in 2–3 × 10−2 sec. Characteristics of the four sample pore spaces, Si , are computed using the methods of Hyman et al. [22] and are reported in Table I. The porosities fall at the lower end of what is considered a normal range, 0.2 n 0.7, and the mean hydraulic radius, tortuosities, and permeabilities correspond to well-sorted gravel or sand [22]. The measurements of tortuosity reported in Table I indicate that fluid particle trajectories deviate from a straight line path, for which τ = 1.
The scalar field of finite time Lyapunov exponents (FTLE) is computed everywhere in the pore space according to the techniques of Shadden et al. [18]. The FTLE, σxT0 , at a given point, x0 , measures the degree of stretching experienced by an infinitesimal fluid element advected with the flow for a finite time T ; larger values indicate more stretching. Specifically, σxT0 = where
1 ln λmax (S(x0 ,T )), 2|T |
S(x0 ,T ) =
dφ T (x(t)) dx(t)
∗
dφ T (x(t)) dx(t)
(2) (3)
is the (right) Cauchy-Green strain tensor, λmax is the maximum eigenvalue of S(x0 ,T ), and A∗ denotes the adjoint of the matrix A. The absolute value is used in (2) to allow for both positive and negative values of T . The strain tensor associated with x0 and T is constructed using the gradient of the flow map φ, which is determined by numerically integrating the trajectory equation, x˙ (x0 ,t) = v(x(x0 ,t)),
(4)
with a fourth order Runge-Kutta scheme coupled with trilinear interpolation to compute the velocities away from computational nodes. The sample pore volume S2 is used to determine a suitable value of T so that σxT0 values have stabilized while still capturing local dynamics, T = 2 × 10−2 sec. The forward flow map, T > 0, is used to identify expansion, while the backward flow map, T < 0, is used to identify contraction. The left side of Fig. 1 is a contour plot of the forward (expanding) FTLE field and the right side is a contour plot of the backward (contracting) FTLE field from one-fourth of a horizontal cross section of sample pore volume S3 (cf. Table I). Contour plots from the other samples are similar. High values in the forward FTLE field form disjoint clusters, and similar characteristics are observed in the backward FTLE field, but those clusters are sparser than clusters observed in the forward FTLE field. This fragmentation indicates that significant expansion and contraction occur in isolated regions within the pore space, and demonstrates the nonuniformity of the flow within pores and across the pore space. To identify locations where significant expansion or contraction in the flow field occur, the forward and backward FTLE fields are thresholded according to σxT0 σ¯ xT + C σˆ xT − σ¯ xT , (5)
063014-2
HYPERBOLIC REGIONS IN FLOWS THROUGH THREE- . . .
PHYSICAL REVIEW E 88, 063014 (2013)
1mm 0
2
4
6
8
10
12
14
16
−3
18 x 10
FIG. 1. (Color online) Contour plots of the forward (left) and backward (right) FTLE fields in one-fourth of a horizontal cross section from a porous medium with porosity 0.38 show that regions of high FTLE values are fragmented. Solid matrix shown in black.
which partitions each field into high and low values [23]. Here, σ¯ xT and σˆ xT are the mean and maximum FTLE values, respectively, and the constant C determines the minimum FTLE value that constitutes significant expansion or contraction. A region of expansion is a connected cluster composed of at least 50 computational grid nodes in the forward FTLE field whose values are greater than the threshold level. Regions of contraction are defined with the same criteria applied to the backward FTLE field. The sample pore volume S2 was used to determine that a value of C = 0.10 extracts the expanding and contracting regions from the respective FTLE fields well. Restricting analysis to clusters composed of at least 50 nodes removes negligible volumes. Table II reports the percentage of the void volume occupied by expanding and contracting regions, the number of regions in each sample, and the correlation coefficient of the expanding and contracting regions with vertical velocity, the primary direction of flow, and the absolute value of horizontal divergence. Expanding regions occupy a larger percentage of the void space than contracting regions indicating that fluid particle trajectories are more prone to disperse than coalesce, and that percentage rises with permeability, porosity, and hydraulic radius. Combined with the decline in the number of regions, this increase indicates that the volume of the expanding regions increases with permeability; images confirm this growth. Across all of the samples, the expanding regions are more
correlated with the vertical velocity and horizontal divergence fields than the contracting regions. In all four samples, the percentage of the void space occupied by the contracting regions is small, less than 20%. The sample with the lowest permeability, S1 , has nearly double the percentage of its void space occupied by these regions compared to the other three samples, indicating that the flow coalesces into select channels through this pore network. This result is consistent with the findings of Hyman et al. [11] who observed that a high percentage of the flux through horizontal cross sections occurred in a low percentage of the void space in low porosity samples. Although the percentage of the void space occupied by the contracting regions remains roughly constant in the higher permeability samples, the number of regions increases. Hence, the regions are smaller in higher permeability samples; images confirm this decrease in volume. The regions of contraction are more correlated with both the vertical velocity and horizontal divergence in lower permeability porous media than in samples with higher permeabilities. The low correlations between the vertical velocity and the absolute value of the horizontal divergence with the expanding and contracting regions suggest that neither of these two quantities of the velocity field alone can fully account for the measured confluence and/or divergence throughout the pore structure. However, because the flow is incompressible, the flow map is volume preserving and the flow field must account for all expansion and contraction. One explanation for these low correlations is that the underlying cause of expansion or contraction depends upon different combinations of the horizontal divergence and vertical velocity and varies spatially. The disparity between locations results from the variable resistance offered by the mesoscopic geometry and topology of the pore structure. The influence of variations in the pore structure on the FTLE fields is demonstrated in Fig. 2 which shows a closeup of the lower right corner of the thresholded forward FTLE field extracted from the left cross section in Fig. 1 at four sequential vertical levels, descending clockwise from the upper left corner. As the height decreases, the pore splits into three and a region of expansion is observed at this junction. Similar behavior is observed in the backward FTLE field where regions of contraction coincide with a convergence in the flow that occurs at pore junctions. IV. HYPERBOLIC REGIONS
TABLE II. Statistics of expanding and contracting regions: sample number (Si ), percent of void volume occupied (PV ), number of regions (NR ), correlation coefficient with vertical velocity (rw ), correlation coefficient with the absolute value of horizontal divergence (rH ). Expansion
Contraction
Si
PV
NR
rw
rH
PV
NR
rw
rH
S1 S2 S3 S4
24.82 27.06 27.76 28.44
176 155 134 137
−0.37 −0.34 −0.31 −0.28
0.36 0.33 0.28 0.23
18.29 8.79 12.24 9.59
251 307 381 389
−0.35 −0.27 −0.28 −0.24
0.33 0.22 0.20 0.12
Hyperbolic regions occur where stable and unstable manifolds in the flow field intersect and are identified by locations that are in both expanding and contracting regions. To provide intuition about the pore geometry responsible for these hyperbolic regions, the left side of Fig. 3 shows a subvolume (about one-thousandth the entire volume) of the pore volume S4 where fluid particle trajectories from two separate pores (red and blue indicate different initial pores) converge into a hyperbolic region (yellow spheres) and then disperse back into the pore network; the primary direction of flow is from top to bottom. The trajectories are separated before they enter the hyperbolic region, intermingle therein, and are mixed together when they disperse into two disjoint pores below.
063014-3
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013) TABLE III. Statistics of hyperbolic regions; sample number (Si ), percent of void volume occupied (PV ), percent of percolating particles that intersect at least one hyperbolic regions (PT ), averages given a particle intersects a hyperbolic region: number of unique intersections (NT ), percent of particle travel time in a hyperbolic region (TT ).
−3
x 10 18 16 14
Unrestricted
Restricted
12
Si
PV
PT
NT
TT
PV
PT
NT
TT
10
S1 S2 S3 S4
9.46 4.81 5.87 4.50
100.00 99.11 100.00 99.07
12.96 10.63 12.06 8.49
17.50 7.07 7.57 5.90
4.43 1.05 1.12 0.61
100.00 65.18 69.49 51.01
3.43 1.88 1.95 1.31
10.64 2.54 2.65 1.98
8 6
this phenomena of in-pore mixing requires that the flow under consideration be three dimensional. In a similar study involving simpler geometry and topology, Cardenas [9] observed vortices near pore throats when considering three dimensional Navier-Stokes flow through spherical and ellipsoidal cubicpacked grains at a range of Reynolds numbers. Table III reports statistics of hyperbolic regions and their intersections with fluid particle trajectories that percolate through the entire domain, determined using the techniques of Hyman et al. [11]. Unrestricted values are generated using all hyperbolic regions, and restricted values are the hyperbolic regions which are composed of at least 50 nodes. The percentages of the void volume occupied by the hyperbolic regions are less than the corresponding values of either the expanding or contracting regions (cf. PV Table II), decrease with increasing permeability, and demonstrate the sparsity of hyperbolic regions, especially those with larger volumes as evidenced by the even smaller percentages reported for the restricted regions. Despite the isolation of the hyperbolic regions that results from the fragmentation of the expanding and contracting regions, they are connected via fluid particle trajectories. Hence, the product of this in-pore mixing is distributed throughout the rest of the pore space via streamlines in the flow. Even though the hyperbolic regions occupy less than
4 2 0
1 mm FIG. 2. (Color online) Closeup of the lower right corner of the thresholded forward FTLE field extracted from the left cross section in Fig. 1 at four sequential vertical levels, descending clockwise from the upper left corner. The thresholded forward FTLE field shows that a region of expansion occurs where a pore splits and the flow separates.
The subfigure on the right is a closeup of the hyperbolic regions on the left with selected trajectories colored yellow. Starting from the black cross in the center of the image, placed at the center of the hyperbolic region, the paths of the yellow trajectories are integrated backwards in time to show where they originate and forward in time to show where they go. The yellow trajectories originate in separate pores above and disperse into two disjoint pores below, thereby demonstrating the pore-scale mixing which occurs in hyperbolic regions. Furthermore, because the flow is laminar and at steady state
1 mm FIG. 3. (Color online) Left: Subvolume of the sample S4 (about one-thousandth the entire volume) with fluid particle trajectories from two separate pores (red and blue indicate different initial pores) converging into a hyperbolic region (yellow spheres) and diverging back into the pore network. Right: Closeup view of the of hyperbolic region on the left. Starting from the black cross in the center of the image, placed at the center of the hyperbolic region, the paths of the yellow trajectories are integrated backwards in time to show where they originate and forward in time to show where they exit the subvolume. 063014-4
HYPERBOLIC REGIONS IN FLOWS THROUGH THREE- . . .
PHYSICAL REVIEW E 88, 063014 (2013)
10% of the void volume, greater than 99% of percolating fluid particles pass through a control volume whose corresponding grid node is in a hyperbolic region. Those that do not remain close to the impermeable lateral boundaries of the domain throughout their trajectories. The lower values reported for the restricted regions indicate that intersections with larger hyperbolic regions are less frequent. Nonetheless, at least half of the particles intersect between one and two of these larger hyperbolic regions, and in the lowest permeability sample all percolating particles still pass through at least one of the larger hyperbolic regions. Moreover, particles that intersect a hyperbolic region pass through ten unique regions on average and spend a nontrivial percent of their travel time therein. Additional statistics not given show that particles intersect a hyperbolic region shortly after insertion. The probability that a particle’s first intersection with a hyperbolic region occurs at a given height decays exponentially with increasing vertical distance from the particle’s insertion point. The first-passage time distributions reported by Hyman et al. [11] illustrate that flows through virtual pore structures sampled from the same ensemble as those used here demonstrate longitudinal dispersion whose magnitude decreases with porosity. The measurements of tortuosity reported in Table I indicate that fluid particle trajectories deviate from a straight line path and indicate that transverse dispersion is also occurring. Computation of particles’ lateral displacement from their initial position shows that on average, particles disperse horizontally at least one-tenth of the domain width as they pass from the top of the domain to the bottom in all of the samples considered. In conjunction with the high percentage of particle intersections with hyperbolic regions (cf. PT in Table III) this suggests that the effects of localized mixing are dispersed laterally as well as in the primary direction of flow through an entire pore structure.
flows through realistic three dimensional pore structures are sparse and disjoint. They coincide with pore throats where both confluence and divergence in the flow occur, and occupy a low percentage of the void volume. Nonetheless, nearly every percolating fluid particle trajectory passes through several hyperbolic regions and together their trajectories form pathways between these mixing regions. The percolating trajectories and hyperbolic regions provide a skeleton of mixing in the flow field and indicate that the effects of localized mixing are distributed across trajectories and throughout an entire pore space. Figure 2 shows a region of expansion occurring where a pore splits, and highlights how the resistance offered by the intricately connected pore structure induces a nonuniform flow therein. Through Lagrangian analysis it was determined that this nonuniformity occurs on two length scales in the flow field. There is a larger scale of nonuniformity across the entire flow field that corresponds to what is classically refereed to as mechanical dispersion. This is evidenced by the contour plots of the forward and backward FTLE fields (Fig. 1), tortuosities of percolating particles (Table I), and their lateral displacement (Sec. IV). There is also a relatively small scale that occurs in individual pores, evidenced by the mixing observed in localized hyperbolic regions (Fig. 3) that has not been extensively investigated previously due to limitations in the ability to observe flows within pore spaces on the one hand, and computational limitations on the other. Advective tracer particles that percolate through the entire domain link the two lengths scales together (Table III) indicating that mixing and transport results from the interaction of these two scales of dispersion.
V. SUMMARY AND REMARKS
ACKNOWLEDGMENTS
Through analysis of finite time Lyapunov exponent fields it was determined that hyperbolic regions in steady state laminar
We would like to thank R. Indik, K. Lin, and M. Tabor for discussions and insightful comments.
[1] J. Bear, Dynamics of Fluids in Porous Media (Dover, New York, 1988). [2] M. Sahimi, Flow and Transport in Porous Media and Fractured Rock (Wiley, Weinheim, 2012). [3] H. Hasimoto, J. Fluid Mech. 5, 317 (1959). [4] B. E. Launder and T. H. Massey, J. Heat Transf. 100, 565 (1978). [5] A. S. Sangani and A. Acrivos, Int. J. Multiphase Flow 8, 343 (1982). [6] A. S. Sangani and A. Acrivos, Int. J. Multiphase Flow 8, 193 (1982). [7] A. Zick and G. M. Homsy, J. Fluid Mech. 115, 13 (1982). [8] A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 54, 406 (1996). [9] M. B. Cardenas, Geophys. Res. Lett. 35, L18402 (2008).
[10] B. Bijeljic, P. Mostaghimi, and M. J. Blunt, Phys. Rev. Lett. 107, 204502 (2011). [11] J. D. Hyman, P. K. Smolarkiewicz, and C. L. Winter, Phys. Rev. E 86, 056701 (2012). [12] B. Bijeljic, A. Raeini, P. Mostaghimi, and M. J. Blunt, Phys. Rev. E 87, 013011 (2013). [13] M. Matyka, Z. Koza, J. Gołembiewski, M. Kostur, and M. Januszewski, Phys. Rev. E 88, 023018 (2013). [14] S. Neuman, Adv. Water Resour. 28, 149 (2005). [15] A. M. Tartakovsky and S. P. Neuman, Geophys. Res. Lett. 35, L21401 (2008). [16] J. M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport (Cambridge University Press, Cambridge, England, 1989), Vol. 3. [17] G. A. Voth, G. Haller, and J. P. Gollub, Phys. Rev. Lett. 88, 254501 (2002).
063014-5
JEFFREY D. HYMAN AND C. LARRABEE WINTER
PHYSICAL REVIEW E 88, 063014 (2013)
[18] S. C. Shadden, F. Lekien, and J. E. Marsden, Physica D 212, 271 (2005). [19] P. K. Smolarkiewicz and C. L. Winter, J. Comput. Phys. 229, 3121 (2010). [20] J. M. Prusa, P. K. Smolarkiewicz, and A. A. Wyszogrodzki, Comput. Fluids 37, 1193 (2008).
[21] R. Mittal and G. Iaccarino, Annu. Rev. Fluid Mech. 37, 239 (2005). [22] J. D. Hyman, P. K. Smolarkiewicz, and C. L. Winter, Water Resour. Res. 49, 2080 (2013). [23] Y. Zhu and R. Newell, Mon. Weather Rev. 126, 725 (1998).
063014-6
