Published in Vadose Zone Journal 3:164-169 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: UNDERSTANDING SUBSURFACE FLOW AND TRANSPORT PROCESSES AT THE IDAHO NATIONAL ENGINEERING & ENVIRONMENTAL LABORATORY (INEEL) SITE |
Particle-Based Direct Numerical Simulation of Contaminant Transport and Deposition in Porous Flow
Ray A. Berry*,
Richard C. Martineau and
Thomas R. Wood
Idaho National Engineering and Environmental Laboratory, P.O. Box 1625, Idaho Falls, ID 83415-2210
Correspondence: * Corresponding author (berry{at}ida.net).
Received for publication 14 May 2003.
 |
ABSTRACT
|
|---|
This work describes an approach to porous flow modeling in which the "micro-length scale to macro-length scale" physical descriptions are addressed as Lagrangian, pore-level flow and transport. The flow features of the physical domain are solved by direct numerical simulation (DNS) with a grid-free, hybrid smoothed particle hydrodynamics (SPH) numerical method (Berry, 2002) based on a local Riemann solution. In addition to being able to handle the large deformation, fluid–fluid and fluid–solid interactions within the contorted geometries of intra- and inter-pore-scale modeling, this Riemann–SPH method should be able to simulate other complexities, such as multiple fluid phases and chemical, particulate, and microbial transport with volumetric and surface reactions. A simple model is presented for the transfer of a contaminant from a carrier fluid to solid surfaces and is demonstrated for flow in a simulated porous media.
Abbreviations: DNS, direct numerical simulations • SPH, smoothed particle hydrodynamics
|
INTRODUCTION
|
|---|
CURRENT TEXTBOOK LITERATURE describes the phenomenological parameters and spatial averaging of material properties and velocity used to generate basic "mixture" continuum field theories for describing and simulating flow in a porous media (e.g., Slattery, 1978; Whitaker, 1999). The many successes of this approach are well documented, but relatively little discussion is devoted to those cases for which predictions with this approach fail. A growing body of evidence is now showing substantial deficiencies in the existing models and predictive capabilities for unsaturated porous flow and transport, particularly in fractured rock vadose zones (Pruess et al., 1999; Glass et al., 2002). In attempts to better describe the transport of water and contaminants through the vadose zone, specialized mixture continuum field theories have been recently extended to include preferential flow, "where water and solutes move along certain pathways, while bypassing a fraction of the porous matrix" (Hendrickx and Flury, 2001), and nonequilibrium flow, where "for various reasons, infiltrating water does not have sufficient time to equilibrate with slowly moving resident water in the bulk of the soil matrix" (Jarvis, 1998). Hendrickx and Flury (2001) discussed preferential flow mechanisms and processes at various scales, ranging from pore to areal scales, while Simunek et al. (2003) gave a review of models describing macropore and nonequilibrium flows resulting from processes at the pore and pedon scales.
The physical description of the disparate length and time scales inherent in geological multiphase fluid transport is traditionally accomplished using continuum theory. This theory views the fluids and solids as interpenetrating mixtures, each being governed by conservation laws either postulated or derived from the method of volume averaging. This Eulerian continuum approach results in field equations for the flow properties of all phases of the multiphase mixture system. As currently practiced, there are at least two major drawbacks with this approach. First, the averaging process produces a loss of information. When manipulating the averaged equations, covariance terms are typically neglected because pore- and fracture-scale models are needed for their construction. The resulting equation system is not complete because the pore- and fracture-scale information is filtered in the averaging process. If these covariance terms are to be included directly, suitable constitutive relations must be developed, in their stead, to obtain closure. Second, averaging leads to unknown terms representing the interactions between the phases (including the solid phase). These terms must be modeled to close the description of the system. The detailed nature of the interactions between solid and fluid phases cannot be understood by the application of mixture theory alone. It is the pore- and fracture-scale physics and the phasic interactions that determine the nature of these constitutive relations. This important information is seldom used, however, because of the difficulty in obtaining it.
As a first step toward an increased understanding of the pore- and fracture-scale physical behavior and phasic interactions, we describe a DNS approach with which a Lagrangian description of fluid flows can be simulated through a broad range of representative porous and/or fracture configurations. The DNS approach provides the capability of studying the nonlinear and geometrically complicated phenomena of multiphase inter- and intra-pore flow wherein preferential and nonequilibrium flow and transport are accounted for as an integral part of the physical description. We present a novel approach that uses a grid-free numerical method to simulate a simplified Lagrangian, microscopic flow. Our approach, however, should extend to include multiple fluid phases with surface tension, liquid bridging, capillary attachment–detachment, chemical constituents with volumetric and surface reactions and pore-surface adsorption–desorption and transport, suspended particulate transport and deposition, and microbial transport, deposition, and buildup. Here we modify the SPH grid-free method, which employs no permanent connectivity between computational nodes or particles, to incorporate local Riemann solutions, to produce our Riemann–SPH grid-free method. It should apply to problems where traditional Eulerian and Lagrangian techniques fail, such as those involving large deformation, fluid–fluid or fluid–solid interactions, and contorted geometries, precisely the situation present in detailed, intra- and inter-pore-scale (and fracture-scale) modeling.
|
RIEMANN–SPH METHOD
|
|---|
The SPH method is a Lagrangian, particle based approach to hydrodynamics. Unlike finite-difference, finite-element, and finite-volume methods, which utilize grid cells to cover the spatial domain of interest, particle methods use computational points which are moving in space where computational data in the material is sampled. The motion of the particles is usually governed by the equations of the material. The SPH method is Lagrangian in the hydrodynamics sense in that the particles follow the material (fluid and perhaps microbes and solids) flow. There is no need for a computational grid because the particles carry all of the computational information, and the particles themselves form the computational framework on which the material flow equations are solved. Originally introduced in the context of astrophysics (Lucy, 1977; Gingold and Monaghan, 1977), the SPH method has evolved into multiple variations with a plethora of applications ranging from high-strain-rate solid deformation mechanics to ice-jammed river flows. This approach offers a substantial advantage for problems that contain moving discontinuities and for those characterized by large deformations as would be encountered in a DNS simulation of saturated and unsaturated, inter- and intra-pore flows. The Lagrangian feature of SPH may also prove ideal for chemical species and biological organism transport, and microscopic surface transport, as well as chemical reaction and biofilm interaction and growth.
With SPH, a moving continuum field is simulated by a large group of moving computational particles. We then replace the point mechanical and material thermomechanical variables (e.g., fluid velocity or fluid pressure at a specific point), which vary rapidly on a scale comparable with the particle spacing, by smoothed variables obtained by averaging over regions that are large compared with the particle spacing but small compared with the complete system. Intuitively, a continuum is modeled as a system of particles where, in the limit, the particles' sizes vanish and the distances between adjacent particles approach zero. In this limit, the system of particles approaches a deformable continuum, that is, an integral over a continuous body that may be considered as the limiting extension of a summation over a corresponding system of particles. A formal mathematical definition of local mean variables is used to translate the point equations for the material flow into computationally attractive representations of the conservation–balance equations.
The basis for SPH is a regularization or smoothing of any point property of the material, f
. The local mean value 
is defined by
 | [1] |
The weighting or smoothing kernel function w
is usually a centrally peaked, spherically symmetric function, with compact support radius h (the smoothing length) representing the effective width of the kernel. The properties of the weighting function are further specified to ensure that the ordinary equations of vector calculus can be performed on the local mean variables, and also that certain other mathematical manipulations, such as differentiating under an integral, can be used. Typical kernels have a Gaussian-like shape, continuous derivatives of specified order, and compact support so that only the volume within a constant multiple of h contributes to the integration. The SPH method is usually viewed as an approximation method for information at unordered points, but it can also be loosely interpreted as a Monte Carlo technique for evaluating Eq. [1] [in which w is a probability distribution function]. A Monte Carlo summation approximation of Eq. [1] is
 | [2] |
where the factor mj/
j is the volume element that SPH associates with particle j, and n is the number of particles within the influence range of the kernel. The spatial gradient can be directly calculated without a computational grid by appropriate differentiation of the approximation formula, Eq. [2], producing
 | [3] |
Divergence of a vector function is found similarly.
The basics of the SPH method are well described in Monaghan (1992). For ease of exposition, the Euler equations with a barotropic equation of state (pressure dependence only on density) will be used to illustrate development of the SPH method on realistic, fluid-dynamic field equations and to show the hybrid extension to incorporate a Riemann solution method. For the porous flow demonstration that follows, viscosity terms are added to produce the Navier–Stokes equations governing viscous Newtonian fluid dynamic motions. A variety of SPH approximations of the inviscid, isothermal hydrodynamic equations (Euler equations) have appeared in the literature. For example,
 | [4] |
and, for the two-dimensional analysis to be demonstrated here, the weighting kernel (Wingate and Fisher, 1993)
is selected. These equations represent a system of ordinary differential equations for each particle that are time-integrated numerically with the explicit-Euler method. Time step size is adjusted dynamically to satisfy the minimum Courant condition, 
t 
0.4
, where c is the sound speed.
Because the basic SPH method exhibits what has become known in the literature as "tensile instability" (can also occur when pressure is positive) (Swegle et al., 1995), and because these compressible formulations also require an artificial viscosity for stabilization, a hybrid SPH method was developed which incorporates a Riemann solution technique. Following the work of Parshikov et al. (2000), for high-velocity impact of solids, we insert (into SPH) approximations for the velocities and pressures obtained from Riemann solutions to break up discontinuities between interacting particles (Berry, 2002). No artificial viscosity is explicitly added, and tensile instability appears to be dramatically reduced or eliminated. In this hybrid approach, the interaction of two particles i and j is taken to be equivalent to that of the acoustic Riemann solution that results from assuming an initial discontinuity between i and j. The mass continuity and momentum balance Eq. [4] become
 | [5] |
and u*Rij and p*ij are normal velocities and pressures resulting from the Riemann solution between particles i and j. Furthermore, the solution between particles i and j is allowed to mimic fluid cavitation or fracture behavior and thereby produce dual values of velocity corresponding to a specified cavitation or fracture pressure, or tension cut-off value (space limitations preclude further discussion here) (Stanyukovich, 1960).
The SPH method described is a compressible method. An approach to represent nearly incompressible flow is given by Rogers et al. (1987). We pick, for example, a stiff, barotropic equation of state such as
 | [6] |
where 0 is the reference density and p0 is the magnitude of the pressure. If these parameters are properly picked to yield (in the simulations) density variations of 0.25 to 1% with corresponding Mach numbers of around 0.1 and a Courant number of about 0.1, then the flow can be regarded as approximately incompressible. This however, may not always be an adequate way to produce the incompressible constraint. For such cases, the method should be modified to update pressure semi-implicitly in a manner consistent with other incompressible fluid dynamic algorithms.
|
CONTAMINANT MODEL
|
|---|
A contaminant transport and interaction model is developed here to show the rich dynamical behavior that can result from even the simplest of physical assumptions. It is assumed that a contaminant may be present in a carrier liquid (water) as it flows, and that this action occurs solely by advection (i.e., without diffusion). All solid surfaces are assumed to have a much higher affinity for the contaminant than does the carrier liquid. Furthermore, they have the ability to pull the contaminant out of the carrier liquid if the liquid is within a prescribed distance, from solid surfaces (thus making contaminant pick-up a localized phenomena). This model is simplified by assuming that, when contaminated fluid is within this prescribed distance of a solid surface, all of that carrier liquid's contaminant is transferred to the solid on a time scale that is very small compared with the transport time scale, and that the carrier liquid is then contaminant-free or clean. For the demonstration example that follows, no limits are placed on how much contaminant may be accumulated by a solid surface and it is assumed that the solid's affinity for the contaminant is so high that it cannot transfer back to the carrier liquid (one-way transfer). It is assumed that contaminated carrier fluid has a contaminant concentration value of 1.0 and that the cleaned carrier fluid has a contaminant concentration of 0.0. Thus, solid surfaces continually accumulate contaminant in unit increments as carrier liquid particles approach solid particles (surfaces) within the prescribed distance for contaminant transference.
|
DEMONSTRATION EXAMPLE
|
|---|
We illustrate this Riemann–SPH method with an example from Berry (2002). Utilizing a particle neighbor-tracking method, edge detection algorithm (Berry, 2001), and basic particle code structures (again space limitations preclude further details here), it is relatively easy to add increasingly more complex physical behaviors. We have constructed a basic model for unsaturated liquid flow through a high-porosity media with a small particulate contaminant transport and deposition model for mock subsurface flows. Though easily extendable to three dimensions, as a first attempt at DNS of porous flows with the Riemann–SPH method, this initial effort solves the Navier–Stokes equations in two dimensions with only one fluid phase and no surface tension. However, it clearly shows the merit of this approach. Because we are ultimately interested in an increased understanding of the pore- and fracture-scale physical behaviors we pick a representative-scale problem to demonstrate.
This simulation models a column of water, initially 1.2 mm deep with hydrostatic pressure head and contaminated with normalized concentration level of 1.0, placed above a column of representative porous material, also 1.2 mm high by 1.0 mm wide but uncontaminated with a normalized value of 0.0. The water is allowed to flow through the porous column under the influence of gravity, and out the bottom, where it falls into a catch-pan and accumulates. The water is represented with 12000 particles, while the solids of the column walls, catch-pan, and the porous region are represented with 5089 particles. The sticking distance for the contaminant to transfer to the solid surfaces is taken to be 0.015 mm. The degree of flow saturation is not only time dependent, but is spatially varying as well. Figure 1
shows the spatial variation of velocity magnitude within the porous column at approximate times of 0.117, 0.230, and 0.459 s. Figure 2
shows the contaminant distribution at the same respective times as those of Fig. 1. Mechanical dispersion effects can be readily observed in these Figures. From Fig. 2c, which shows the distribution of contaminant when the porous column has nearly drained, it is apparent that a significant portion of the contaminant has been sorbed onto the porous solid particles, thereby reducing its presence in the pan's water. It is noted for Fig. 2 that the color scale of the contaminant level indicates normalized values of 0.0 (uncontaminated, shown as blue), 1.0 (contaminated, shown black), and 3.0 (contaminated with normalized values of 3.0 or higher, shown as magenta), while the contaminant level on some of the porous solid particles approaches normalized values of 150. (Normalized values of 3 or 150 indicate contaminant concentrations 3 times or 150 times, respectively, higher than the initial contaminations concentration level in the water column.) This demonstration problem was executed for approximately 5.5 h on a PC Athlon 1.3-GHz computer, completing roughly 40000 time steps. Though the time step sizes were dynamically adjusted, typical time steps were on the order of 12 µs.
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 2. Contaminant distribution at same three successive times as Fig. 1: (a) 0.117 s, (b) 0.230 s, and (c) 0.459 s. Black indicates normalized contaminate level of 1.0, magenta a normalized contaminant level of 3.0 or higher. Blue indicates uncontaminated fluid or solid (normalized contaminant level of 0.0).
|
|
A simulation of this type would be very difficult or impossible using traditional numerical methods. To apply an Eulerian approach, a fixed, boundary-fitted grid would first have to be set up, sophisticated boundary conditions applied, and a very good interface-tracking algorithm incorporated. In addition to setting up the initial, boundary-fitted grid, because of the contorted flow paths, a Lagrangian grid-based method would require extensive remeshing, which is computationally expensive and introduces significant error. The simple, history-based, contaminant transport–deposition model exhibited here would be problematic for an Eulerian, grid-based solution method. Whenever fluid or contaminant behavior depends on the history of what that fluid or contaminant has seen in its dynamical evolution, Eulerian methods are at a significant disadvantage.
|
CONCLUSIONS
|
|---|
The simple demonstration example presented here shows the benefits of pursuing the Lagrangian, direct numerical simulation approach with the Riemann–SPH method. We have presented an approach for simulating the physics and chemistry of microscopic flow in porous and/or fractured media. Furthermore, we advocate this approach as a means to building the macroscopic averaged flow models and the appropriate constitutive relations necessary to provide the closure for these macroscopic models. A simple, dynamical contaminant model was shown to produce a result rich in structure that would be difficult to obtain with conventional methods. To extend its applicability to the complexities of multiphase flow, we recommend modification of the method to treat pressure implicitly (in time), addition of another set of particles to represent the gaseous phase, and the incorporation surface tension. We further recommend improvement of the models to incorporate diffusion and chemical reactions, and potentially solid matrix compressibility, deformation, and erosion.
|
ACKNOWLEDGMENTS
|
|---|
The authors would like to acknowledge the U.S. Department of Energy, Office of Environmental Management, for funding this research as part of the Environmental Systems Research Program during Fiscal-Year 2001. This research was conducted under DOE Idaho Operations Office Contract DE-AC07-99ID13727.
|
REFERENCES
|
|---|
- Berry, R.A. 2001. Environmental systems research FY 2000 annual report. p. 218. In D.L. Miller et al. (ed.) Idaho National Engineering and Environmental Laboratory Report INEEL/EXT-01-01085. Available online at http://www.inel.gov/esra/reports.shtml (verified 25 Nov. 2003). INEEL, Idaho Falls, ID.
- Berry, R.A. 2002. Environmental systems research FY 2001 annual report. p. 157. In D.L. Miller et al. (ed.) Idaho National Engineering and Environmental Laboratory Report INEEL/EXT-01-01402. Available online at http://www.inel.gov/esra/reports.shtml (verified 25 Nov. 2003). INEEL, Idaho Falls, ID.
- Gingold, R.A., and J.J. Monaghan. 1977. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society 181:375.
- Glass, R.J., M.J. Nicholl, A.L. Ramirez, and W.D. Daily. 2002. Liquid phase structure within an unsaturated fracture network beneath a surface infiltration event: Field experiment. Water Resour. Res. 38:1199.
- Hendrickx, J.M.H., and M. Flury. 2001. Uniform and preferential flow, mechanisms in the vadose zone, conceptual models of flow and transport in the fractured vadose zone. National Research Council, National Academy Press, Washington, DC.
- Jarvis, N.J. 1998. Modeling the impact of preferential flow on nonpoint source pollution. p. 195–221. In H.M. Selim and L. Ma (ed.) Physical nonequilibrium in soils: Modeling and application. Ann Arbor Press, Chelsea, MI.
- Lucy, L.B. 1977. A numerical approach to the testing of the fission hypothesis. Astron. J. 82:1013.
- Monaghan, J.J. 1992. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys. 30:543.
- Parshikov, A.N., et al. 2000. Improvements in SPH method by means of interparticle contact algorithm and analysis of perforation tests at moderate projectile velocities. Int. J. Impact Eng. 24:779.
- Pruess, K., B. Faybishenko, and G.S. Bodvarsson. 1999. Alternative concepts and approaches for modeling flow and transport in thick unsaturated zones of fractured rocks. J. Contam. Hydrol. 38:281–322.
- Rogers, S.E., D. Kwak, and U. Kaul. 1987. On the accuracy of the pseudocompressibility method in solving the Navier-Stokes equations. Appl. Math. Model. 11:35.
- Simunek, J., et al. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272:14–35.
- Slattery, J.C. 1978. Momentum, Energy, and Mass Transfer in Continua. Krieger Publ. Co., New York.
- Stanyukovich, K.P. 1960. Unsteady motion of continuous media. Pergamon Press, London.
- Swegle, J., D.L. Hicks, and J. Attaway. 1995. Smoothed particle hydrodynamics stability analysis. J. Comp. Phys. 116:123.
- Whitaker, S. 1999. The method of volume averaging. Kluwer Academic Publishers, Dordrecht, The Netherlands.
- Wingate, C.A., and H.N. Fisher. 1993. Strength modeling in SPHC. Los Alamos National Laboratory Report LA-UR-93-3942. LANL, Los Alamos, NM.
Copyright © 2009 by Soil Science Society of America
TOP
ABSTRACT
INTRODUCTION
