Published online 16 August 2005
Published in Vadose Zone J 4:848-855 (2005)
DOI: 10.2136/vzj2004.0178
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Simulation of Unsaturated Flow in Complex Fractures Using Smoothed Particle Hydrodynamics
Alexandre M. Tartakovskya,* and
Paul Meakinb
a Pacific Northwest National Lab., P.O. Box 999/MS K1-85, Richland, WA 99352
b Idaho National Lab., P.O. Box 1625, MS 2025, Idaho Falls, ID 83415-2025
* Corresponding author (Alexandre.Tartakovsky{at}pnl.gov)
Received 9 December 2004.
 |
ABSTRACT
|
|---|
Smoothed particle hydrodynamics (SPH) models were used to simulate unsaturated flow in fractures with complex geometries. The SPH is a fully Lagrangian particle-based method that allows the dynamics of interfaces separating fluids to be modeled without employing complex front tracking schemes. In SPH simulations, the fluid density field is represented by a superposition of weighting functions centered on particles which represent the fluids. The pressure is related to the fluid density through an equation of state, and the particles move in response to the pressure gradient. The SPH does not require the construction of grids that would otherwise introduce numerical dispersion. The model can be used to simulate complex free-surface flow phenomenon such as invasion of wetting and nonwetting fluids into three-dimensional fractures. These processes are a severe challenge for grid-based methods. Surface tension was simulated by using a van der Waals equation of state and a combination of short-range repulsive and longer-range attractive interactions between fluid particles. The wetting behavior was simulated using similar interactions between mobile fluid particles and stationary boundary particles. The fracture geometry was generated from self-affine fractal surfaces. The fractal model was based on a large body of experimental work, which indicates that fracture surfaces have a self-affine fractal geometry characterized by a material independent (quasi universal) Hurst exponent of about 0.75.
Abbreviations: SPH, smoothed particle hydrodynamics
|
INTRODUCTION
|
|---|
UNSATURATED FLOW IN FRACTURES plays an important role in the migration of fluids in the subsurface under natural and altered conditions. For example, fractures may significantly decrease the time required for contaminants to migrate through the vadose zone to an underlying aquifer. Computer modeling is playing an ever increasing role in the development of a better understanding of fluids and the prediction of fluid flow in oil reservoirs, hydrothermal reservoirs, contaminated sites, aquifers, and other systems. However, the application of standard grid-based numerical methods to free-surface and multiphase fluid flow processes with complex dynamical interfaces is fraught with difficulties such as artificial interface broadening and grid entanglement. An alternative approach is to use particles to represent the fluids. The fluid–air or fluid–fluid interfaces move with the particles, there is no need to explicitly track the interfaces, and processes such as fluid fragmentation and coalescence can be handled without difficulty.
In this paper we describe a numerical model based on SPH to simulate three-dimensional invasion of nonwetting and wetting fluids in fracture apertures with complex geometries. Such complex processes are very difficult to model with traditional grid-based numerical schemes. The SPH is an interpolation-based technique that can be used to solve numerically systems of partial differential equations. The Lagrangian particle nature of SPH allows physical and chemical effects to be incorporated into the modeling of flow processes, which may also be complicated by irregular and deformable boundaries. The simplicity of the SPH method allows complex flow problems to be solved with relative ease in one, two, and three dimensions. The SPH was first introduced by Lucy (1977) and Gingold and Monaghan (1977) for the solution of Navier–Stokes equations in the context of astrophysical fluid dynamics. Since its introduction, SPH has been successfully used to model a wide range of fluid flow processes and the behavior of solids subjected to large deformations. For example, Monaghan (1994) used SPH to model the collapse of dams, Morris et al. (1997) extended SPH to model low Reynolds number flows, and Zhu et al. (1999) and Zhu and Fox (2001) applied SPH to study pore-scale flow and transport in saturated porous media. Incorporation of the effects of surface tension into SPH simulations has been a vexing problem. Morris (2000) modeled surface tension based on its continuous macroscopic description with surface tension forces that are proportional to the fluid–fluid interface curvature. This approach gives an accurate estimation of the effects of surface tension but involves rather complex calculations of front curvatures that in some cases may lead to significant errors. Nugent and Posch (2000) used attractive forces, corresponding to the cohesive pressure in the van der Waals equation of state to simulate surface tension in two-dimensional SPH simulations. They found that it was necessary to increase the range of the attractive forces to at least twice the range of the SPH weighting function to obtain stable vapor drops, a costly solution causing a significant increase in a computational time. In the work described in this paper we used a different approach. Initially we used van der Waals equation of state to simulate formation and behavior of three-dimensional fluid drops. When we used a low particle density this resulted in the formation of a liquid drops that did not have a stable smooth spherical shape, a problem similar to that encountered by Nugent and Posch (2000). When we increased the particle density (that corresponds to increase in resolution of simulation) we were able to create stable spherical fluid drops but it also led to collapse of some particles. When a combination of short-range repulsive and large-range attractive interactions (with the range of the attractive interactions equal to the range of the SPH weighting function) were added to the three-dimensional SPH equations spherical liquid drops were formed without collapsing particles, and this approach allows surface tension, phase separation, and different wetting conditions to be modeled. The use of a combination of short-range repulsive and long-range attractive interactions was motivated by the molecular origins of surface tension.
|
SMOOTHED PARTICLE HYDRODYNAMICS REPRESENTATION OF FLUID FLOW EQUATIONS
|
|---|
Smoothed particle hydrodynamics theory has been very well described by Monaghan (1992). Here we present a short summary. Any continuous field, A(r), at position r can be approximated by
 | [1] |
where the integration is performed over the entire space of the field and the weighting function W with a support scale of 6h (a range of 3h) becomes a Dirac delta function in the h 
0 limit (Fig. 1)
. The weighting function satisfies the normalization condition
 | [2] |
In SPH, the fluid is represented by a discrete set of N particles. The position of the ith particles is denoted by the vector ri, i = 1,...,N, and each particle has a mass, mi, and density, 
i. The properties associated with any particle i are calculated by approximating the integral in Eq. [1] by the sum
 | [3] |
and the magnitude of the field at position r is approximated by
 | [4] |
where the summation is performed over all of the particles.
A variety of forms, including spline functions of different order, have been used for the weighting functions. The order of the spline function can be viewed as being analogous to the order of finite element discretization schemes. We found that for accurate representation of free surfaces a spline function of at least fourth order is needed, and we used the fourth-order weighting function
 | [5] |
where 
= 
and = 
in two and three spatial dimensions.
Here, we present an application of SPH to free-surface flows in fracture apertures with complex geometries. Such flows are governed by the Navier–Stokes equations describing the conservation of momentum and conservation of mass:
 | [6] |
and
 | [7] |
where v is the continuous fluid velocity vector, P and are the continuous fluid pressure and density fields, 
is the dynamic viscosity, and g is the gravitational acceleration vector.
In SPH Eq. [6] and [7] become equations of motion for each particle, and fluid flow is represented by the movement of the particles. In the simplest versions of SPH, all of the particles have the same mass, and conservation of the number of particles rigorously conserves the mass of the fluid.
A number of SPH approximations to the Navier–Stokes equation have been described in the literature (Monaghan, 1992). We used the momentum conservation equation
 | [8] |
defined for each particle. In this equation, the SPH representation of the pressure gradient (used in the first term on the right hand site of Eq. [8]) was derived by Monaghan (1992) and the SPH representation of the viscous force was obtained by Morris et al. (1997).
Equation [7] determines the evolution of the particle density. For simplicity, and exact mass conservation, the density of each particle is commonly computed directly from Eq. [3] with A being the density (Nugent and Posch, 2000; Monaghan, 1992):
 | [9] |
Given the densities a, the pressure in Eq. [8] is obtained from the equation of state. In this work we used a van der Waals equations of state in which the pressure is given by (Nugent and Posch, 2000)
 | [10] |
where 
= kb/m (kb is the Boltzmann constant), 
1 = c1/m and 2 = c2/m. Here, c1 and c2 are the van der Waals constants, and m is the mass of the particles.
To simulate the effect of surface tension we used a van der Waals equation of state and added the effects of particle–particle interactions into Eq. [8] with interaction forces given by
 | [11] |
where sij is the strength of the force acting between particles i and j (Fig. 2)
. The total force due to interparticle interactions acting on any particle i then can be found from
 | [12] |
Since Fij = –Fji, the particle–particle interactions conserve momentum and can be safely added into Eq. [8]. The range of the interaction force coincides with the range of the weighting function for maximum computational efficiency. For a given fluid, the magnitude of this force depends only on the distance between particles. The force is repulsive for distances less than h, attractive for distances between h and 3h and is zero for distances larger than 3h. Apart from the effects of small density and configurational fluctuations in the interior of the fluids, the total attractive force acting on the particles is nonzero only at the fluid–fluid interface and next to the walls of the fracture. In the boundary region, this force acts in the direction of the density gradient creating surface tension.
We verified the accuracy of the surface tension model by studying large gravity free oscillations of three-dimensional liquid drops that were initially deformed into oblate (discus-like) and prolate (cigar-like) shapes. In the van der Waals equation of state, the dimensionless values of the constants used were: kT = 0.067, 1 = 0.014 and 2 = 0.0052. The strength of interaction force was sij = 0.01. The dimensionless viscosity was set to zero and dimensionless masses of the particles were set to unity. To prepare the initial liquid drops, we first placed particles on a regular cubical lattice inside a sphere with radius equals 5 in units of the range, 3h, of the weighting function. An SPH simulation was then run to allow the particle system to reach equilibrium in the presence of the free surface. The resulting drop was then slowly deformed between two reflecting plates into an oblate shape or it was squeezed inside a reflecting cylinder with slowly decreasing radius to form a prolate shape. When the constraints of the partially confining reflective surfaces were released, the drops underwent oscillations depicted in Fig. 3 and 4
. The oscillations of the initially oblate drop closely resemble the oscillations of a golf ball size liquid drop observed experimentally in the space shuttle Columbia (Apfel et al., 1997). Even though dynamic viscosity is zero, the amplitudes of oscillations are decreasing due to the intrinsic viscosity that is inherent to particle systems (Posch et al., 1995).
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 3. Gravity-free oscillation of a fluid sphere initially deformed into an oblate shape. First row is a side view and second row is a top view.
|
|
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 4. Gravity-free oscillations of a fluid sphere initially deformed into a prolate shape. First row is a side view and second row is a top view.
|
|
Another advantage of using particle–particle interactions Fij is the ability to model a variety of wetting and nonwetting behaviors at the solid boundaries. This can be achieved by placing stationary particles (denoted by subscript b) in the vicinity of the boundaries and by assigning different interaction strengths swb and snb between the nonwetting, n, and wetting, w, fluid particles and the boundary (fracture wall) particles. The strength of the wetting is determined by the strength of the fluid-boundary particle interactions relative to the strength of the fluid–fluid interactions. Because of the short-range repulsive part of Fij, the fluid–boundary interparticle interactions also simulate no-flow boundary condition in a manner similar to that employed by Monaghan (1994), who used the repulsive part of a Lennard–Jones interaction acting between the boundary and fluids particles. Occasionally fluid particles will overcome the repulsive interactions with the wall particles. These particles are prevented from penetrating the fracture walls by using bounce-back boundary conditions in combination with the interaction with boundary particles. Figure 5
shows the transient velocity profile for two-dimensional flow between two parallel plates. The solid lines are velocity profiles obtained analytically, and asterisks and circles represent velocity profiles obtained from SPH simulations with two different approaches to modeling the solid boundary. The van der Waals Eq. [10] was used in these examples to close the SPH equations with dimensionless values of the constants: kT = 0.08, 1 = 0.017 and 2 = 0.0075. The dimensionless viscosity was 8 and the dimensionless masses of the particles were set to unity. In these simulations the fluid density was high enough to prevent phase separation. The asterisks show the results of an SPH simulation using bounce back boundary condition in combination with short-range repulsive and long-range attractive particle–particle interaction forces. The strength of the interaction force was sij = 0.001. The circles indicate the results obtained using the approach proposed by Morris et al. (1997) for treating no-slip boundary conditions at solid surfaces. This approach consists of assigning artificial velocities vb given by
 | [13] |
to the boundary particles b for each fluid particle i in the vicinity of the solid boundary.
In Eq. [13], db and di are the normal (shortest) distances from the solid boundary to the boundary particle and the fluid particle. The artificial boundary particle velocities are then included in the calculation of the viscous forces in Eq. [8], and this procedure accurately simulates impermeable no-slip conditions at solid boundaries. Even though the bounce back method of treating solid boundaries leads to velocity fields that deviate more from analytical results than the velocity fields obtained with the artificial velocity method, the simplicity of implementation and computational efficiency justify the usage of bounce back boundary conditions in more complex models for our present purposes.
A linked-list approach with an underlying cubic lattice (with the size of the lattice unit equal to the range of the weighting function, 3h, for computational efficiency) was used to rapidly locate all of the particles within a range 3h of any selected particle.
Flow in Fracture Apertures
To simulate flow in three-dimensional fractures the aperture geometry was generated from self-affine fractal surfaces. The fractal model is based on extensive observations that the fracture surfaces of brittle materials, including rocks (Schmittbuhl et al., 1995), have a self-affine fractal geometry (Mandelbrot, 1984; Mandelbrot et al., 1984), which can be characterized by a Hurst exponent with a more or less material independent (quasi universal) value of about 0.75 (Bouchaud, 1997; Bouchaud et al., 1990). The fracture aperture was simulated as the gap between a self-affine surface with a Hurst exponent of 0.7 and a replica of the surface, which was translated both horizontally and vertically, without rotation and with periodic boundary boundaries in the directions parallel to the plane of the fracture. Figure 6
shows the resulting fracture aperture field that ranges from 0.5 to 9 in the units of the range, 3h, of the weighting function. The fracture size in the lateral directions (parallel to the plane of the fracture) was 32 x 32 in units of the range, 3 h, of the weighting function (Fig. 7)
and the (vertical) thickness of the computational domain was 16.
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 7. Sketch of a three-dimensional fracture. Fluid is injected along one edge of the fracture (at x = 0). Gravity acts in the x direction.
|
|
Interactions with the boundaries of the fracture were modeled using ‘boundary’ particles that are immobile but were included in the calculation of the total force acting on fluid particles in Eq. [8] and [12]. Initially particles were placed randomly in the 32 by 32 by 16 box and then given time to relax. Once particles reached a uniform distribution, the fracture geometry was imposed on the particle system. The particles within a unit distance (3h) of self-affine surfaces were ‘frozen’ and labeled as boundary particles creating the fracture walls. The rest of the particles were removed, except for particles with x coordinate less than unity. These particles created a ‘fluid pool’ at the entrance of fracture. The flow was initiated by applying a gravitational force. A constant head boundary condition was applied at the entrance of the fracture (x = 0) by keeping a constant number of particle in the pool. At each time step, particles were added to or (if necessary) removed from the fluid pool at randomly selected positions to compensate exactly for the net flow of particles into the aperture. To prevent particles from leaving the system in the direction opposite to flow, bounce back boundary conditions were used at x = 0. At each time step, the densities at each of the particles were calculated using Eq. [9] and the pressure at each particle was obtained using the equation of state [10]. The acceleration of each fluid particle, ai = dvi/dt, was calculated from Eq. [8] with the forces calculated from Eq. [11] and [12] added to the total force in Eq. [8]. The new particle position was then calculated using the "velocity Verlet" algorithm (Allen and Tildesley, 2001), which takes the form
 | [14] |
and
 | [15] |
In the van der Waals equation of state, the dimensionless values of the constants used were: kT = 0.1, 1 = 0.022 and 2 = 0.013. The dimensionless viscosity and dimensionless masses of the particles were set to unity.
Free Surface Flow of a Nonwetting Fluid in a Fracture Aperture
Figures 8 and 9
show the invasion of a nonwetting fluid into an empty fracture aperture for different gravity fields (g = 0.01 and g = 0.0001, respectively, where g is the component of acceleration due to gravity acting in the x direction). The X-Z cross-section (cross-sections perpendicular to the plane of the fracture in the direction of flow) as well as a view looking down onto the X-Y plane are shown at three stages during the simulation. The wetting conditions were modeled by making snn = 0.01 for the interactions between the particles representing the nonwetting fluid and snb = 0.001 for the interactions between the fluid particles and the boundary particles. These figures show that the nonwetting fluid forms a convex front in the direction of the flow, as would be expected for a nonwetting fluid. The convex shape of the front is more pronounced for large gravitational fields g than for smaller g. For a smaller gravitational field our model predicts that the fracture will be completely filled with nonwetting fluid while for a larger gravitational field the fracture becomes only partially filled.
View larger version (36K):
[in this window]
[in a new window]
|
Fig. 8. Injection of a nonwetting (light particles) fluid into a fracture. g = 0.01. The X-Z cross-section and the X-Y top view are shown at three different times. The dark particles are immobile particles modeling the walls of the fracture.
|
|
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 9. Injection of a nonwetting (light particles) fluid into a fracture. g = 0.0001. The X-Z cross-section and the X-Y top view are shown at three different times. The dark particles are immobile particles modeling the walls of the fracture.
|
|
Free Surface Flow of a Wetting Fluid in the Fracture
The invasion of a wetting fluid into a fracture was modeled by making the interaction strength sww between wetting fluid particles equal to 0.01 and the interaction strength swb between fluid and boundary particles equal to 0.02. Figures 10 and 11
show the flow for g = 0.01 and 0.001. These figures show that a concave front is formed in the direction of the flow. The concave shape is more pronounced for a flow under the influence of a smaller gravitational field. This phenomenon (an increasing of contact angle with increasing velocity of advancing front) has been also observed experimentally (for example by Schwartz and Tejeda, 1972).
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 10. Injection of a wetting (light particles) fluid into a fracture. g = 0.01. The X-Z cross-section and the X-Y top view are shown at three different times. The dark particles are immobile particles modeling the walls of the fracture.
|
|
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 11. Injection of a wetting (light particles) fluid into a fracture. g = 0.001. The X-Z cross-section and the X-Y top view are shown at three different times. The dark particles are immobile particles modeling the walls of the fracture.
|
|
|
SUMMARY
|
|---|
The use of particle–particle interactions in combination with traditional SPH methods allows complex fluid behavior, including the large oscillation of fluid drops, the effects of surface tension and different fluid wetting behaviors, to be realistically simulated. The complex dynamic of free surfaces, which would be very difficult to simulate with grid-based methods, can be simulated without undue effort. The power of Lagrangian particle-based methods such as SPH is illustrated by these simulations. Detailed comparisons with laboratory experiments and results obtained from well-established numerical methods are needed to quantitatively validate the SPH-based modeling methods described in this paper.
|
ACKNOWLEDGMENTS
|
|---|
This work was supported by the Laboratory Directed Research and Development program at the Idaho National Laboratory (INL). The first author was partially supported by the Laboratory Directed Research and Development program at the Pacific Northwest National Laboratory (PNNL). The INL is operated for the USDOE by Battelle under Contract DE-AC07-05ID14517. The PNNL is operated for the USDOE by Battelle under Contract DE-AC06-76RL01830.
|
REFERENCES
|
|---|
- Apfel, R.E., Y. Tian, J. Jankovsky, T. Shi, X. Chen, R.G. Holt, E. Trinh, A. Croonquist, K. Thornton, A. Sacco, Jr., C. Coleman, R. Leslie, and D. Matthiesen. 1997. Free oscillations and surfactant studies of superdeformed drops in microgravity. Phys. Rev. Lett. 78:1912–1915.[CrossRef]
- Allen, M.P., and D.J. Tildesley. 2001. Computer simulation of liquids. Oxford Univ. Press, Oxford, UK.
- Bouchaud, E. 1997. Scaling properties of cracks. J. Physics-Condensed Matter 9:4319–4344.[CrossRef]
- Bouchaud, E., G. Lapasset, and J. Planes. 1990. Fractal dimension of fractured surfaces—A universal value. Europhys. Lett. 13:73–79.
- Gingold, R.A., and J.J. Monaghan. 1977. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181:375–389.
- Lucy, L.B. 1977. Numerical approach to the testing of the fission hypothesis. Astron. J. 82:1013–1024.[CrossRef]
- Mandelbrot, B.B. 1984. The fractal geometry of nature. W.H. Freeman, New York.
- Mandelbrot, B.B., D.E. Passoja, and A.J. Paullay. 1984. Fractal character of fracture surfaces of metals. Nature (London) 308:721–722.[CrossRef]
- Monaghan, J.J. 1994. Simulating free surface flows with SPH. J. Comp. Phys. 110:399–406.[CrossRef]
- Monaghan, J.J. 1992. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys. 30:543–574.[CrossRef]
- Morris, J.P., P.J. Fox, and Y. Zhu. 1997. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136:214–226.[CrossRef]
- Morris, J.P. 2000. Simulating surface tension with smoothed particle hydrodynamics. Int. J. Numer. Methods Fluids 33:333–353.[CrossRef]
- Nugent, S., and H.A. Posch. 2000. Liquid drops and surface tension with smoothed particle applied mechanics. Phys. Rev. E 62:4968–4975.[CrossRef]
- Posch, H.A., W.G. Hoover, and O. Kum. 1995. Steady-state shear flows via nonequilibrium molecular-dynamics and smooth-particle applied mechanics. Phys. Rev. E 52:1711–1720.[CrossRef]
- Schmittbuhl, J., F. Schmitt, and C. Scholz. 1995. Scaling invariance of crack surfaces. J.Geophys. Res.-Solid Earth 100:5953–5973.
- Schwartz, A.M., and S.B. Tejeda. 1972. Studies of dynamic contact angles on solids. J. Colloid Interface Sci. 38:359–375.[CrossRef]
- Zhu, Y., P.J. Fox, and J.P. Morris. 1999. A pore-scale numerical model for flow through porous media. Int. J. Numer. Anal. Methods Geomechanics 23:881–904.[CrossRef]
- Zhu, Y., and P.J. Fox. 2001. Smoothed particle hydrodynamics model for diffusion through porous media. Transp. Porous Media 43:441–471.[CrossRef]
TOP
ABSTRACT
INTRODUCTION
