(210a) Inertial Effects on Dispersion in Porous Media

            In 1905 Charles Slichter made both field and laboratory measurements of dispersion in porous media, and correctly identified the variation of the velocity field as the primary cause of dispersion in porous media [1].  However, it was not until Taylor's seminal paper on dispersion in a capillary tube in 1953 [2] that the understanding of hydrodynamic dispersion was put into a sound theoretical context.  Since that time, a wide array of theoretical approaches has been used for predicting (i.e., upscaling) the effective pore-scale dispersion tensor in porous media.  In the research reported here, we have focused on the influence of inertial contributions to the dispersion tensor in porous media at moderate Reynolds numbers.

            Although it has been recognized for some time that inertial effects can influence the observed conductivity for fluid flow through porous media, the influence that inertial flows have on dispersion has not been as widely realized.  In industrial processes, the flow through a packed bed reactor frequently has a particle Reynolds number well within the inertial range [3].  Specific examples of such processes include gas adsorption, catalysis, filtration, bed washing operations, heat transfer, combustion, environmental hydraulics, and even immobilized enzyme reactors.

            For the work reported here, we have adopted the method of volume averaging [4] with closure for predicting the effective dispersion tensor from the microscale geometry and fluid properties.  Volume averaging upscales heterogeneous media by forming spatial averages of the balance equations of the independent variables of interest.  To start out with, on formulates the conservation equations of interest at the sub-pore scale.  In our case, we assume incompressible flow, and the sub-pore scale mass and momentum balance equations take the form (where we will follow the method and notation described by Whitaker [4])

Here, c_A&gamma is the concentration of chemical species A in the &gamma-phase (the fluid phase), v_&gamma is the microscale velocity vector of the fluid phase, D_A&gamma is the diffusion coefficient for a dilute solution of species A in the fluid, &rho is the fluid density, m is the fluid viscosity, and p_&gamma is the fluid pressure.  The averaging volume is indicated by V, and the fluid-filled volume is indicated by V_&gamma.  The volume fraction of fluid is given by &epsilon_&gamma = V/V_&gamma.  The notation A_&gamma&kappa indicates the fluid-solid interface, and n_&gamma&kappa is the unit normal to this surface directed outward from the fluid phase.

            The expressions given by the microscale equations above can be upscaled using the superficial and intrinsic volume averages are denoted by (respectively)

            Following thethe developments described by Whitaker [4], one can apply these averages to develop a macroscale mass balance equation of the form

the effective dispersion tensor takes the form

where where the vector b_&gamma is interpreted as the integral of a Greens function of the 'closure problem' (see Appendix B in ref. [5]).  The closure problem is set of ancillary equations that are solved over a representative region of the porous medium.  The purpose of the closure problem is to determine the effective macroscale influence of the microscale physics and structure; extensive discussion of the closure problem is given elsewhere (6).

            For the closure of the problem posed here, we have adopted the simple unit cell shown in Fig 1(a).  The closure problem was solved numerically to determine the effective longitudinal dispersion coefficient. The commercial finite-element code COMSOL Multiphysics was used to determine the solution; convergence was determined using a heuristic convergence analysis via grid refinement.  The parameters used for these simulations were as follows: &epsilon_&gamma =0.376, &rho = 1000 kg/m³, and D_A&gamma = 1x10^-9 m²/s.  The length of one side of the unit cell was set at l = 1x10^-3 m.  Only a quarter of the domain needed to be simulated because of symmetry, and for this region approximately 10,000 elements were used to resolve the solutions.  In the solution to the equations for the concentration deviations, a streamline upwind Petrov-Galerkin method was used to reduce oscillations.

            Results from these simulations are given in Fig. 2.  In this figure, the normalized effective longitudinal dispersion coefficient is plotted versus the particle Peclet and Reynolds numbers

In Fig. 2, it is apparent that the results from the simple unit cell appearing in Fig. 1 generally underpredicts the effective longitudinally hydrodynamic dispersion coefficient for a fixed value of the Peclet number.  The underprediction observed is most likely the result of (1) the simple structure of the unit cell (i.e., no longer-range structures involving multiple particles are represented, and these would generally increase dispersion), and (2) numerical dispersion, which would make the actual Peclet number smaller than intended (resulting in a shift of the curve to the right).  Importantly, however, the slope of the curve is correct, including the concave downward shape that is observed for the data.  For Peclet numbers above about 1, the relationship between the effective longitudinal dispersion coefficient and the Peclet number is given by a power law of the form

For the results found in this study, the exponent is in the range 0.95 < &beta < 1.2, with &beta generally decreasing with increasing Peclet number.  This result consistent with the experimental data.

Figure 1.  Structure of the flow field for two Reynolds numbers for a periodic array.  (a) The section of the periodic array shown.  (b) Stokes flow for Re_p = 0.1.  (c) Inertial effects are evident at moderate Reynolds numbers (Re_p = 200), where both jet and vortex formation can be seen.

Figure 2. The effective longitudinal dispersion coefficient as a function of the Peclet number: Comparison of closure simulations and experimental data.  The inclusion of inertial effects provides one possible explanation for the change in slope of this curve at high values of the Peclet (or, equivalently, Reynolds) number.

References

[1]       Slichter, C. S. Field measurement of the rate of movement of underground waters, US Geological Survey, 1905.

[2]       Taylor, S. G. Proc. Roy. Soc. London 1953, 219, 186-203.

[3]       Strang, D. A.; Geankoplis, C. J. Ind. Eng. Chem. 1958, 50, 1305-1308.

[4]       Whitaker, S. The Method of Volume Averaging; Kluwer: Dordrecht, 1999.

[5]       Golfier, F.; Quintard, M.; Cherblanc, F.; Harvey, C. F.; Zinn, B.; Wood, B. D. Submitted to Advances in Water Resources 2006.