- Email: [email protected]
Lattice Boltzmann simulations and contact angle hysteresis in convergent–divergent media
Journal of Petroleum Science and Engineering 33 (2002) 161 – 171 www.elsevier.com/locate/jpetscieng
Lattice Boltzmann simulations and contact angle hysteresis in convergent–divergent media R.D. Hazlett*, R.N. Vaidya Potential Research Solutions, 1818 Shelmire Drive, Dallas, TX 75224, USA
Abstract The lattice Boltzmann method was used to perform multiphase fluid dynamics calculations that include preferential fluid interactions with solid walls. Contact angle was calibrated with respect to wall interaction parameters using captive drops and observed LaPlace pressure differences across interfaces. In simulations, captive drops were forced to propagate down straight tubes to observe dynamic contact angles. Captive drops were also propagated down sinusoidally constricted capillaries in simulations to observe dynamic contact angles in convergent – divergent systems. Rate effects on observed contact angle were also probed. Negligible hysteresis or rate effect was observed over the conditions examined, indicating that the phenomenon of contact angle hysteresis is not a result of hydrodynamic interactions with an equilibrium wetting condition. It is conjectured that the time scale for solid – fluid interaction is too short in the current lattice Boltzmann scheme in comparison with that for hydrodynamic interaction. The results suggest that such lattice Boltzmann models must incorporate dynamic rock – fluid interactions different from static ones in order to portray these known phenomena or slow the mathematical equilibration process for preferential wetting. D 2002 Published by Elsevier Science B.V. Keywords: Boltzmann; Contact angle; Hysteresis; Simulation
1. Introduction 1.1. Lattice Boltzmann method 1.1.1. General overview The lattice Boltzmann method (LBM) is a computational fluid dynamics technique which tracks the populations of fluid particles on the continuum scale as they interact with each other and rigid walls, whilst conserving mass and momentum. The underlying assumption of LBM is that fluid flow behavior is not sensitive to fluid microstructure, and thus, a fluid dy*
Corresponding author. E-mail address: [email protected] (R.D. Hazlett).
namic method can be constructed on the continuum level with packets of interacting fluid. Single-phase interactions follow an ideal gas model, and Navier – Stokes behavior is recovered in the limit of low Mach number (Wolfram, 1986; Frisch et al., 1986). The method is well established. The LBM has been extended to multiphase flow with viscosity contrast and interfacial phenomena considerations (Rothman and Keller, 1988; Gunstensen et al., 1991; Grunau et al., 1993). Readers should refer to the work of these authors for detail on the technique or comparison with other methods. In contrast to three-phase line motion studies using molecular dynamics (Collet et al., 1997), this work investigates the ability of the LBM to model such hysteretic wetting angle phenomena.
0920-4105/02/$ - see front matter D 2002 Published by Elsevier Science B.V. PII: S 0 9 2 0 - 4 1 0 5 ( 0 1 ) 0 0 1 8 3 - 8
162
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
1.1.2. Interfacial tension emulation While other methods have been researched (Orlandini et al., 1995; Langaas, 1999), the mode of interfacial tension and wettability emulation used in these investigations is an adaptation of the Gunstensen scheme (Gunstensen et al., 1991). This scheme is observed to produce thin, stable interfaces and recovers the LaPlace pressure jump across interfaces with the introduction of an interfacial tension strength parameter, A. The interfacial tension in the LBM is the slope of a plot between Boltzmann pressure jump and inverse drop size. The method invokes a threestage process (Grunau, 1993). In the first stage, interfacial cells are identified as those containing significant populations of more than one fluid color or type. Secondly, the preferred interfacial orientation is identified for those cells. Population densities are examined locally. A planar orientation is selected through a minimization process which would require the least amount of work to segregate fluid particles of different color. Once chosen, the second stage, known as recoloring, commences. Diffusing particle populations in the local neighborhood are returned toward their parent phase to the extent controlled by the interfacial tension, A. This is illustrated in Fig. 1a and is repeated every time step. 1.1.3. Wettability emulation Wettability, herein defined as the competitive interaction of two or more fluid phases with a solid, is a special case of interfacial tension emulation. Immobile fluid particles are placed within the outer shell of solids with populations of phases in proportion to desired wetting characteristics. As embedded particle populations, the stationary mass acts as attractors and strongly influences the process of preferential interfacial orientation at the three-phase line of contact. This is illustrated in Fig. 1b. An important point to note is that interfacial positioning, and thus contact angle, is determined from a local computational procedure involving information only from nearest and next nearest neighbor cells in a cubic lattice. While the Young – Dupre´ relationship, and variants thereof (Good, 1952; Hazlett, 1990), can be derived from thermodynamics involving total contact areas, global minimization methods would not be appropriate for a LBM. In LBM, the apparent contact angle
Fig. 1. Interfacial tension and wettability modeling within the lattice Boltzmann method as adapted from Gunstensen: (a) IFT and phase separation result from mass redistribution with respect to a locally preferred interfacial orientation, and (b) stationary mass is imposed on the solid boundary, which influences interfacial orientation near walls.
results from local enforcement of the LBM rules for fluid – fluid –solid interaction using an interfacial tension parameter value which satisfies the macroscopic LaPlace equation (Grunau, 1993). The method of
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
surface interaction bias through the insertion of colored rest mass is somewhat akin to alteration of the surface free energy of the solid relative to contacting fluid. It could also be envisioned as defining the presence of a permanent wetting film, although no hydrodynamics are portrayed within this layer of ill-defined thickness. Since each surface site can acquire coexisting color populations, the full range of wetting behavior and wettability character can be achieved. This procedure for wettability emulation, however, requires calibration. 1.2. Contact angle hysteresis 1.2.1. Evidence While the equilibrium contact angle has a firm thermodynamic basis, it is actually difficult to observe because of mechanical equilibrium considerations. Experiments follow procedures which introduce time history as a variable. Propagation of a three-phase line of contact is inevitable, and contact angle hysteresis is observed. Good (1952) ascribed the phenomenon of contact angle hysteresis to the existence of an energy barrier which represents ‘‘the contortional energy that the drop must have if the drop is to move from one favored configuration of the triple interface (more or less meandering across the ‘hills and valleys’ of the solid surface) to another.’’ A good illustration of the contortions which must take place for interfacial ad-
163
vancement is presented by Oliver and Mason (1977). Fig. 2 shows photomicrographs of fluid advancement on a grid with heterogeneity comparable to drop size. Good was obviously referring to the influence of surface roughness, but hysteresis can similarly result from chemically heterogeneous surfaces (Hazlett, 1992). There appears to be a close analogy between relationships developed for surface roughness and chemical heterogeneity (de Gennes, 1985). Sticking or pinning (Joanny and de Gennes, 1984) of the interface can occur at surface defects, which causes entire ranges of apparent contact angles to be observed. Thus, the equilibrium contact angle is observed only in the most carefully controlled conditions. A further distinction between advancing and receding contact angles is made in an attempt to categorize the three-phase line motion history. As sample purity, inhomogeneity, equilibration time, and interfacial advancement rate also influence observed contact angles, the reporting of only advancing and receding angles should be accompanied with supplementary information. A velocity dependence of contact angle has been reported by numerous researchers (Hoffman, 1975; Cain et al., 1983). When considering propagating drops or ganglia, fluid and three-phase line motion at leading and trailing edges become coupled events (Payatakes and Dias, 1984). The existence of contact angle hysteresis is indisputable; the reason for the observed hysteresis is not always obvious.
Fig. 2. A look at the origin of contact angle hysteresis and metastable states during three-phase line movement after Oliver and Mason (1977): (a) fluid on a very coarse heterogeneous gridded surface, (b) the next low energy state during drop expansion from the previous photo, (c) closeup examination of the three-phase line of contact indicating considerable distortion in drop shape.
164
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
1.2.2. Controversy Two main controversies exist in the literature on the mechanisms involved in modeling three-phase line motion, and thus contact angle hysteresis. The initial point of contention is the global versus local enforcement of physical laws. Is Young– Dupre´ to be enforced locally? That is, does the equilibrium contact angle prevail at the three-phase line and observed behavior result from global surface minimization, or must models employ some other physical law or rule set on the microscopic level. A second issue involves the hydrodynamic boundary conditions that most accurately reflect contact line dynamics. Line motion has been modeled as both slippage across the surface and as interfacial rolling. The rolling hypothesis would involve interfacial motion, much like that of tank track. Oliver and Mason’s (1977) experiments resembled three-phase line slippage. Other works have invoked slippage at the boundary to match or model motions (Dussan, 1976; Miksis and Davis, 1994). The dye experiments of Dussan and Davis (1974), however, were more supportive of the interfacial rolling mechanism. The LBM employed in this work takes a stand on both issues. As described previously, the LBM exclusively uses local interactions in the bulk and at interfaces. Thus, local enforcement of Young – Dupre´ controls observed behavior. This enforcement, however, is at the continuum level, and not microscopic. Still, surface roughness and inhomogeneity are described at this same level of detail, so that surface defects influence local interfacial orientation computations. Although the LBM identifies interfaces and performs additional computations to preserve interfacial tension and integrity, it does not employ a separate suite of hydrodynamic interaction rules at such sites. The LBM used also does not allow exchange of mass between fluid and solid surface cells. Thus, LBM modeling of three-phase line motion resembles three-phase line slippage rather than rolling. Fluid – solid hydrodynamic interaction is modeled in LBM with a first-order accurate bounce-back condition analogous to no-slip for shear stress on the wall (Grunau, 1993). This is not to be confused with a no-slip boundary condition for the interface, since the interface position is determined each time step following all hydrodynamic interactions with the walls. Higher order accuracy in LBM solid interaction is still an active topic of research. Results
are presented within the context of these underlying modeling assumptions.
2. Methodology 2.1. Boundary conditions 2.1.1. Geometry This study of contact angle and hysteresis effects made use of simple geometry in contrast to application-oriented studies performed by the authors (Hazlett et al., 1998; Vaidya et al., 1996) in imaged reservoir rock pore networks. Two types of virtual 3-D conduit geometry were used in this work: (a) straight capillary tubes of varying diameter and (b) sinusoidally constricted straight capillaries. Subroutines were added to the LBM code to generate tube geometry with userdefined minimum diameter, maximum diameter, and wavelength. Straight capillary tubes were a degenerate case. The ability for the LBM to accurately portray hydrodynamics in these geometries is supported by Fig. 3, which shows a comparison of the simulated and theoretical permeability for each geometry with selected parameters. Theoretical values were taken from the work of Neira and Payatakes (1979). 2.1.2. Wettability calibration In the calibration of wettability, a captive drop method was developed. A straight capillary tube was generated. The desired wettability was applied by setting wall color and density. A slug of one phase was centrally placed within the capillary with fluid of opposite type on either side. The density of all phases was uniformly set. No pressure gradient was imposed on exterior ends of the capillary. The fluid –fluid – solid system was allowed to evolve to an equilibrium condition. While visual inspection was available for contact angle interpretation, a more accurate method was used. A central longitudinal density profile was examined for the capillary pressure jump at each interface. The pressure difference, along with knowledge of the capillary internal diameter, allowed the determination of cosh through the Young – Dupre´ equation. An example profile is provided in Fig. 4. Tests were repeated for different wall colors, densities, and drop color in order to construct a calibra-
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
165
Fig. 6, indicating a size effect exists for capillaries under 5 lattice units in diameter, and that it persists in capillaries below 9 lattice units in diameter for wettability assignments near the sensitive region of the calibration curve. 2.1.3. External boundary conditions The LBM rules for fluid– fluid –solid interaction are set once fluid parameters are prescribed, solid walls are identified, and wetting conditions applied. Still, users must set external boundary conditions to initiate organized flow. In essence, the LBM is merely a suite of interaction rules that must be followed once the system is stimulated externally. Thus, we must additionally define initial and external boundary conditions. Many LBM users like to invoke forcing functions. These represent body forces and may be appropriate for gravity-driven flows. This work made use of an externally controlled pressure driving force in fluid propagation studies. In the LBM, this is achieved by setting a lattice Boltzmann density gradient. 2.2. Numerical experiments 2.2.1. Spontaneous imbibition In order to probe for process sensitivity to applied wetting conditions, spontaneous imbibition experiments were performed. Without gravity present, capillaries were filled with nonwetting fluid of prescribed
Fig. 3. Comparison of theoretical and simulated permeability for straight and sinusoidally constricted capillaries: (a) permeability and aspect ratio for straight capillaries as a function of diameter, and (b) theoretical versus simulated permeability for a series of sinusoidally constricted capillaries of the same maximum diameter and wavelength, but differing minimum constriction diameter.
tion for wall density versus contact angle. The developed calibration is provided in Fig. 5 for an 11cell capillary diameter. Note the extreme sensitivity of the contact angle to the wall color density difference near F 1.8. This calibration was tested for smaller capillary sizes. The results are furnished in
Fig. 4. Examining the central longitudinal pressure profile to extract the contact angle from the LaPlace pressure jump at each fluid – fluid interface. The graphical inset shows a central cross-section of the phase map within the 3-D tube.
166
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
Fig. 5. Contact angle calibration as determined by the captive drop technique as a function of the assigned wetting density difference. Graphical inserts depict a central slice of the phase distribution within the capillary for the various points. Note that the contact angle reported is that measured through the blue (gray in these images) phase.
uniform density. At time zero, a wetting phase was allowed to contact the capillary end. Co-current imbibition of wetting phase was monitored with time.
2.2.2. Hysteresis experiments In hysteresis studies, captive drops were forced to propagate under the influence of a defined external density gradient. Wettability parameters were chosen based upon the equilibrium contact angle calibration. For propagating drops, the same pressure jump technique was used although a viscous pressure drop was superimposed. Advancing and receding angles were analyzed. Tests were conducted for different density gradients to probe for velocity effects. Similar tests were performed in sinusoidally constricted capillaries.
3. Results and discussion 3.1. Spontaneous imbibition in cylindrical capillaries
Fig. 6. Sensitivity of wettability assignment with respect to capillary size. Contact angles measured through the blue phase.
Spontaneous imbibition was observed to be sensitive to the equilibrium contact angle assigned to the walls. This sensitivity persisted beyond the point of perfect wetting. A plot of spontaneous imbibition into cylindrical capillaries is given in Fig. 7 as a function
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
167
Fig. 7. Sensitivity to wettability as indicated by spontaneous imbibition rate into a cylindrical capillary. Note that sensitivity to the wetting strength still occurs below the point of perfect wetting (h = 0). Contact angles shown are measured through the blue phase.
of wettability assignment. Although not compared directly to theoretical or experimental imbibition rates, the demonstrated sensitivity is in agreement with expectations.
faces does not reflect the same time or length scale for hydrodynamic interactions. 3.3. Propagation in sinusoidally constricted capillaries
3.2. Propagation in cylindrical capillaries Fig. 8 shows a plot of pressure profile for propagating captive drops in a cylindrical capillary for different driving forces, and thus velocities. Mass velocities were computed and used to distinguish curves in the legend. Within statistical accuracy, there was no difference observed between advancing and trailing contact angles, nor was evidence of velocity dependence displayed. These experiments were performed in a capillary of only 5 lattice units in diameter. With the walls in such close proximity, the drop shape must have been strongly influenced by conditions imposed at the walls. Still, the inability for hydrodynamics to affect the interfacial curvature is a signal that perhaps the rigor in which the equilibrium condition is applied to construct and preserve inter-
Further tests were performed for more complex, sinusoidally constricted capillaries. Without straight walls, it was anticipated to be a composite effect of equilibrium contact angle applied locally and the oscillating slope of the solid wall with respect to the tube axis. For fixed maximum diameter, wavelength, and equilibrium contact angle, propagating slugs were monitored for two different choices of minimum constriction diameter. Results are shown in Figs. 9 and 10. Clearly, interfacial motion for advancing and trailing edges becomes correlated, although the results should vary widely with slug size relative to medium wavelength. The apparent contact angles at both advancing and receding interfaces take on a range of values, challenging the meaning of apparent contact angle for convergent –divergent porous media. In both cases,
168
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
Fig. 8. Examining the central longitudinal pressure profile to extract advancing and receding contact angles, with respect to the blue phase, from captive drops propagating down a cylindrical capillary tube. No statistically significant contact angle hysteresis was observed, nor was velocity dependence evidenced over the range studied. Wetting parameter set at 2.5 to yield an equilibrium angle of 37j.
Fig. 9. Examining central longitudinal pressure jumps at leading and trailing interfaces as a function of time for a captive drop propagation down a sinusoidally constricted tube. Capillary parameters are as indicated. Wetting set for he = 65j.
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
169
Fig. 10. Examining central longitudinal pressure jumps at leading and trailing interfaces as a function of time for captive drop propagation down a sinusoidally constricted tube. Minimum capillary diameter reduced to 3 lattice units. Wetting set for he = 65j.
the maximum pressure difference is observed when the equilibrium angle is found macroscopically. Both advancing and receding angles are found to achieve this maximum. In Fig. 9, the trailing interface is even observed to take on negative pressure differences, indicating apparent wetting behavior counter to the naturally imposed wettability. This behavior could be easily predicted for capillary-dominated flow conditions as a hemispherical cap oscillates in curvature due to the combination of local equilibrium contact angle and wall curvature.
4. Conclusions 4.1. General observations The lattice Boltzmann method offers unique testing opportunities with a higher degree of control and analysis capability than physical experiments on wetting behavior. It was demonstrated that the Gunstensen method can be calibrated to yield local contact angle assignments. It also was observed that different
wettability assignments produced different kinetic behavior for numerical spontaneous imbibition experiments. 4.2. Hysteresis in LBM 4.2.1. Observations Hysteresis effects were briefly examined by simulating a propagating captive drop in a virtual capillary tube of defined equilibrium contact angle. Analysis of experiments readily yielded accurate estimates of dynamic contact angle through the use of the pressure profile across the interfaces. In the limited testing performed, no statistically significant differences were observed between advancing and receding contact angles in straight tubes. In addition, no velocity dependence of contact angles was observed. 4.2.2. Recommendations The inability to make hysteresis and velocitydependent observations challenges the current mode of interfacial orientation determination and preservation in the LBM. Although seemingly rigorous in its
170
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
approach, the Gunstensen method performed every time step may not allow enough interplay between interfacial physics and hydrodynamics. While other methods are under investigation, introduction of a relaxation parameter may allow a balancing of the time scale for surface interactions and that for hydrodynamic interactions. Angles other than the equilibrium contact angle could additionally be imposed to force simulated behavior to conform to experimental observation, but adding ad hoc calibrations is undesirable. Rather, efforts should be made to portray the underlying physics from a mathematical perspective and allow the power of a physically meaningful fluid dynamics tool to integrate the physics over time and space. The real challenges are not in conducting flows in simple geometry, but the ones that are so complex that intuition no longer gives target behavior (Hazlett et al., 1998; Vaidya et al., 1996). Additional testing with single interfaces and other capillary sizes or shapes is in order.
t F k LBM Pc
4.3. Flow in sinusoidally constricted capillaries
The authors gratefully acknowledge the support and assistance of Drs. S.Y. Chen, D.W. Grunau, and W.E. Soll in this collaborative effort. A portion of the work was conducted using the resources of the Advanced Computer Laboratory at Los Alamos National Laboratory.
Additional simulations in sinusoidally constricted tubes indicated advancing and trailing interfaces were changing curvature resulting from interplay between a locally imposed equilibrium contact angle and the slope of the walls with respect to the central axis. If wall curvature were neglected, interpretation of the pressure jump across interfaces would be interpreted as variations in apparent contact angle. This type of analysis is not new, but the interaction between capillarity and hydrodynamics to look at hysteresis in a controlled environment with LBM was novel. The observed curvature variation does not lessen the criticism of the wettability portrayal method. It does, however, illustrate that the LBM can be used to study apparent contact angle in surface roughness or heterogeneous wettability environments. Nomenclature A interfacial tension controlling parameter d Boltzmann density, mass/volume dmax maximum capillary diameter, lattice units dmin minimum capillary diameter, lattice units db Boltzmann blue phase density, mass/volume dr Boltzmann red phase density, mass/volume D capillary diameter, lattice units
S Sw ux / k r h he ha hr
surface normal medium permeability, length2 lattice Boltzmann method capillary pressure or LaPlace pressure difference wetting strength parameter, db dr at the solid surface wetting phase saturation, volume/volume Boltzmann mass velocity, lattice units/time step medium porosity, volume/volume periodic medium wavelength, l.u. interfacial tension contact angle equilibrium contact angle advancing contact angle receding contact angle
Acknowledgements
References Cain, J.B., Francis, D.W., Venter, R.D., Neumann, A.W., 1983. Dynamic contact angles on smooth and rough surfaces. J. Colloid Interface Sci. 94 (1), 123 – 130. Collet, P., De Coninck, J., Dunlop, F., Regnard, A., 1997. Dynamics of the contact line: contact angle hysteresis. Phys. Rev. Lett. 79 (19), 3704 – 3707. de Gennes, P.G., 1985. Wetting: statics and dynamics. Rev. Mod. Phys. 57 (3/1), 827 – 863. Dussan, V.E.B., 1976. The moving contact line: the slip boundary condition. J. Fluid Mech. 77, 665 – 684. Dussan, V.E.B., Davis, S., 1974. On the motion of a fluid – fluid interface along a solid surface. J. Fluid Mech. 65, 71 – 95. Frisch, U., Hasslacher, B., Pomeau, Y., 1986. Lattice – gas automata for the Navier – Stokes equation. Phys. Rev. Lett. 56 (14), 1505 – 1508. Good, R.J., 1952. A thermodynamic derivation of Wenzel’s modification of Young’s equation for contact angles; together with a theory of hysteresis. J. Am. Chem. Soc. 74, 5041 – 5042. Grunau, D.W., 1993. Lattice Methods for Modeling of Hydrodynamics. Ph.D. Thesis, Colorado State University.
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171 Grunau, D.W., Chen, S., Eggert, K., 1993. A lattice Boltzmann model for multiphase fluid flows. Phys. Fluids A 5 (10), 2557 – 2562. Gunstensen, A.K., Rothman, D.H., Zaleski, D.H., Zanetti, G., 1991. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43 (8), 4320 – 4327. Hazlett, R.D., 1990. Fractal applications: wettability and contact angle. J. Colloid Interface Sci. 137 (2), 527 – 533. Hazlett, R.D., 1992. On surface roughness effects in wetting phenomena. J. Adhes. Sci. Technol. 6 (6), 625 – 633. Hazlett, R.D., Chen, S.Y., Soll, W.E., 1998. Wettability and rate effects on immiscible displacement: lattice Boltzmann simulation in microtomographic images of reservoir rocks. J. Pet. Sci. Eng. 20 (3/4), 167 – 175. Hoffman, R., 1975. J. Colloid Interface Sci. 50, 228. Joanny, J.F., de Gennes, P.G., 1984. A model for contact angle hysteresis. J. Chem. Phys. 81, 552 – 562. Langaas, K., 1999. Modelling of Immiscible Two-Phase Flow in Porous Media with the Binary Fluid Lattice Boltzmann Method. Ph.D. Thesis, University of Bergen, Norway. Miksis, M.J., Davis, S.H., 1994. Slip over rough and coated surfaces. J. Fluid Mech. 273, 125 – 139.
171
Oliver, J.F., Mason, S.G., 1977. Microscopic studies on rough surfaces by scanning electron microscopy. J. Colloid Interface Sci. 60, 480 – 487. Orlandini, E., Swift, M.R., Yeomans, J.M., 1995. A lattice Boltzmann model of binary-fluid mixtures. Europhys. Lett. 32 (6), 463 – 468. Neira, M.A., Payatakes, A.C., 1979. Collocation solution of creeping Newtonian flow through sinusoidal tubes. AIChE J. 25, 725 – 730. Payatakes, A.C., Dias, M.M., 1984. Immiscible microdisplacement and ganglion dynamics in porous media. Rev. Chem. Eng. 2 (2), 85 – 174. Rothman, D.H., Keller, J.M., 1988. Immiscible cellular-automaton fluids. J. Stat. Phys. 52 (3/4), 1119 – 1127. Vaidya, R.N, Hazlett, R.D., Buckles, J.J., Coles, M.E., Grunau, D.W., Chen, S.Y., Eggert, K.G., 1996. Simulation of wettability effects on displacement processes by cellular automata methods. In: Morrow, N.R. (Ed.), Proceedings 3rd International Symposium on Evaluation of Reservoir Wettability and Its Effect on Oil Recovery. University of Wyoming, Laramie, pp. 89 – 93. Wolfram, S., 1986. Cellular automaton fluids: I. Basic theory. J. Stat. Phys. 45 (3/4), 471 – 526.
Lattice Boltzmann simulations and contact angle hysteresis in convergent–divergent media R.D. Hazlett*, R.N. Vaidya Potential Research Solutions, 1818 Shelmire Drive, Dallas, TX 75224, USA
Abstract The lattice Boltzmann method was used to perform multiphase fluid dynamics calculations that include preferential fluid interactions with solid walls. Contact angle was calibrated with respect to wall interaction parameters using captive drops and observed LaPlace pressure differences across interfaces. In simulations, captive drops were forced to propagate down straight tubes to observe dynamic contact angles. Captive drops were also propagated down sinusoidally constricted capillaries in simulations to observe dynamic contact angles in convergent – divergent systems. Rate effects on observed contact angle were also probed. Negligible hysteresis or rate effect was observed over the conditions examined, indicating that the phenomenon of contact angle hysteresis is not a result of hydrodynamic interactions with an equilibrium wetting condition. It is conjectured that the time scale for solid – fluid interaction is too short in the current lattice Boltzmann scheme in comparison with that for hydrodynamic interaction. The results suggest that such lattice Boltzmann models must incorporate dynamic rock – fluid interactions different from static ones in order to portray these known phenomena or slow the mathematical equilibration process for preferential wetting. D 2002 Published by Elsevier Science B.V. Keywords: Boltzmann; Contact angle; Hysteresis; Simulation
1. Introduction 1.1. Lattice Boltzmann method 1.1.1. General overview The lattice Boltzmann method (LBM) is a computational fluid dynamics technique which tracks the populations of fluid particles on the continuum scale as they interact with each other and rigid walls, whilst conserving mass and momentum. The underlying assumption of LBM is that fluid flow behavior is not sensitive to fluid microstructure, and thus, a fluid dy*
Corresponding author. E-mail address: [email protected] (R.D. Hazlett).
namic method can be constructed on the continuum level with packets of interacting fluid. Single-phase interactions follow an ideal gas model, and Navier – Stokes behavior is recovered in the limit of low Mach number (Wolfram, 1986; Frisch et al., 1986). The method is well established. The LBM has been extended to multiphase flow with viscosity contrast and interfacial phenomena considerations (Rothman and Keller, 1988; Gunstensen et al., 1991; Grunau et al., 1993). Readers should refer to the work of these authors for detail on the technique or comparison with other methods. In contrast to three-phase line motion studies using molecular dynamics (Collet et al., 1997), this work investigates the ability of the LBM to model such hysteretic wetting angle phenomena.
0920-4105/02/$ - see front matter D 2002 Published by Elsevier Science B.V. PII: S 0 9 2 0 - 4 1 0 5 ( 0 1 ) 0 0 1 8 3 - 8
162
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
1.1.2. Interfacial tension emulation While other methods have been researched (Orlandini et al., 1995; Langaas, 1999), the mode of interfacial tension and wettability emulation used in these investigations is an adaptation of the Gunstensen scheme (Gunstensen et al., 1991). This scheme is observed to produce thin, stable interfaces and recovers the LaPlace pressure jump across interfaces with the introduction of an interfacial tension strength parameter, A. The interfacial tension in the LBM is the slope of a plot between Boltzmann pressure jump and inverse drop size. The method invokes a threestage process (Grunau, 1993). In the first stage, interfacial cells are identified as those containing significant populations of more than one fluid color or type. Secondly, the preferred interfacial orientation is identified for those cells. Population densities are examined locally. A planar orientation is selected through a minimization process which would require the least amount of work to segregate fluid particles of different color. Once chosen, the second stage, known as recoloring, commences. Diffusing particle populations in the local neighborhood are returned toward their parent phase to the extent controlled by the interfacial tension, A. This is illustrated in Fig. 1a and is repeated every time step. 1.1.3. Wettability emulation Wettability, herein defined as the competitive interaction of two or more fluid phases with a solid, is a special case of interfacial tension emulation. Immobile fluid particles are placed within the outer shell of solids with populations of phases in proportion to desired wetting characteristics. As embedded particle populations, the stationary mass acts as attractors and strongly influences the process of preferential interfacial orientation at the three-phase line of contact. This is illustrated in Fig. 1b. An important point to note is that interfacial positioning, and thus contact angle, is determined from a local computational procedure involving information only from nearest and next nearest neighbor cells in a cubic lattice. While the Young – Dupre´ relationship, and variants thereof (Good, 1952; Hazlett, 1990), can be derived from thermodynamics involving total contact areas, global minimization methods would not be appropriate for a LBM. In LBM, the apparent contact angle
Fig. 1. Interfacial tension and wettability modeling within the lattice Boltzmann method as adapted from Gunstensen: (a) IFT and phase separation result from mass redistribution with respect to a locally preferred interfacial orientation, and (b) stationary mass is imposed on the solid boundary, which influences interfacial orientation near walls.
results from local enforcement of the LBM rules for fluid – fluid –solid interaction using an interfacial tension parameter value which satisfies the macroscopic LaPlace equation (Grunau, 1993). The method of
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
surface interaction bias through the insertion of colored rest mass is somewhat akin to alteration of the surface free energy of the solid relative to contacting fluid. It could also be envisioned as defining the presence of a permanent wetting film, although no hydrodynamics are portrayed within this layer of ill-defined thickness. Since each surface site can acquire coexisting color populations, the full range of wetting behavior and wettability character can be achieved. This procedure for wettability emulation, however, requires calibration. 1.2. Contact angle hysteresis 1.2.1. Evidence While the equilibrium contact angle has a firm thermodynamic basis, it is actually difficult to observe because of mechanical equilibrium considerations. Experiments follow procedures which introduce time history as a variable. Propagation of a three-phase line of contact is inevitable, and contact angle hysteresis is observed. Good (1952) ascribed the phenomenon of contact angle hysteresis to the existence of an energy barrier which represents ‘‘the contortional energy that the drop must have if the drop is to move from one favored configuration of the triple interface (more or less meandering across the ‘hills and valleys’ of the solid surface) to another.’’ A good illustration of the contortions which must take place for interfacial ad-
163
vancement is presented by Oliver and Mason (1977). Fig. 2 shows photomicrographs of fluid advancement on a grid with heterogeneity comparable to drop size. Good was obviously referring to the influence of surface roughness, but hysteresis can similarly result from chemically heterogeneous surfaces (Hazlett, 1992). There appears to be a close analogy between relationships developed for surface roughness and chemical heterogeneity (de Gennes, 1985). Sticking or pinning (Joanny and de Gennes, 1984) of the interface can occur at surface defects, which causes entire ranges of apparent contact angles to be observed. Thus, the equilibrium contact angle is observed only in the most carefully controlled conditions. A further distinction between advancing and receding contact angles is made in an attempt to categorize the three-phase line motion history. As sample purity, inhomogeneity, equilibration time, and interfacial advancement rate also influence observed contact angles, the reporting of only advancing and receding angles should be accompanied with supplementary information. A velocity dependence of contact angle has been reported by numerous researchers (Hoffman, 1975; Cain et al., 1983). When considering propagating drops or ganglia, fluid and three-phase line motion at leading and trailing edges become coupled events (Payatakes and Dias, 1984). The existence of contact angle hysteresis is indisputable; the reason for the observed hysteresis is not always obvious.
Fig. 2. A look at the origin of contact angle hysteresis and metastable states during three-phase line movement after Oliver and Mason (1977): (a) fluid on a very coarse heterogeneous gridded surface, (b) the next low energy state during drop expansion from the previous photo, (c) closeup examination of the three-phase line of contact indicating considerable distortion in drop shape.
164
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
1.2.2. Controversy Two main controversies exist in the literature on the mechanisms involved in modeling three-phase line motion, and thus contact angle hysteresis. The initial point of contention is the global versus local enforcement of physical laws. Is Young– Dupre´ to be enforced locally? That is, does the equilibrium contact angle prevail at the three-phase line and observed behavior result from global surface minimization, or must models employ some other physical law or rule set on the microscopic level. A second issue involves the hydrodynamic boundary conditions that most accurately reflect contact line dynamics. Line motion has been modeled as both slippage across the surface and as interfacial rolling. The rolling hypothesis would involve interfacial motion, much like that of tank track. Oliver and Mason’s (1977) experiments resembled three-phase line slippage. Other works have invoked slippage at the boundary to match or model motions (Dussan, 1976; Miksis and Davis, 1994). The dye experiments of Dussan and Davis (1974), however, were more supportive of the interfacial rolling mechanism. The LBM employed in this work takes a stand on both issues. As described previously, the LBM exclusively uses local interactions in the bulk and at interfaces. Thus, local enforcement of Young – Dupre´ controls observed behavior. This enforcement, however, is at the continuum level, and not microscopic. Still, surface roughness and inhomogeneity are described at this same level of detail, so that surface defects influence local interfacial orientation computations. Although the LBM identifies interfaces and performs additional computations to preserve interfacial tension and integrity, it does not employ a separate suite of hydrodynamic interaction rules at such sites. The LBM used also does not allow exchange of mass between fluid and solid surface cells. Thus, LBM modeling of three-phase line motion resembles three-phase line slippage rather than rolling. Fluid – solid hydrodynamic interaction is modeled in LBM with a first-order accurate bounce-back condition analogous to no-slip for shear stress on the wall (Grunau, 1993). This is not to be confused with a no-slip boundary condition for the interface, since the interface position is determined each time step following all hydrodynamic interactions with the walls. Higher order accuracy in LBM solid interaction is still an active topic of research. Results
are presented within the context of these underlying modeling assumptions.
2. Methodology 2.1. Boundary conditions 2.1.1. Geometry This study of contact angle and hysteresis effects made use of simple geometry in contrast to application-oriented studies performed by the authors (Hazlett et al., 1998; Vaidya et al., 1996) in imaged reservoir rock pore networks. Two types of virtual 3-D conduit geometry were used in this work: (a) straight capillary tubes of varying diameter and (b) sinusoidally constricted straight capillaries. Subroutines were added to the LBM code to generate tube geometry with userdefined minimum diameter, maximum diameter, and wavelength. Straight capillary tubes were a degenerate case. The ability for the LBM to accurately portray hydrodynamics in these geometries is supported by Fig. 3, which shows a comparison of the simulated and theoretical permeability for each geometry with selected parameters. Theoretical values were taken from the work of Neira and Payatakes (1979). 2.1.2. Wettability calibration In the calibration of wettability, a captive drop method was developed. A straight capillary tube was generated. The desired wettability was applied by setting wall color and density. A slug of one phase was centrally placed within the capillary with fluid of opposite type on either side. The density of all phases was uniformly set. No pressure gradient was imposed on exterior ends of the capillary. The fluid –fluid – solid system was allowed to evolve to an equilibrium condition. While visual inspection was available for contact angle interpretation, a more accurate method was used. A central longitudinal density profile was examined for the capillary pressure jump at each interface. The pressure difference, along with knowledge of the capillary internal diameter, allowed the determination of cosh through the Young – Dupre´ equation. An example profile is provided in Fig. 4. Tests were repeated for different wall colors, densities, and drop color in order to construct a calibra-
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
165
Fig. 6, indicating a size effect exists for capillaries under 5 lattice units in diameter, and that it persists in capillaries below 9 lattice units in diameter for wettability assignments near the sensitive region of the calibration curve. 2.1.3. External boundary conditions The LBM rules for fluid– fluid –solid interaction are set once fluid parameters are prescribed, solid walls are identified, and wetting conditions applied. Still, users must set external boundary conditions to initiate organized flow. In essence, the LBM is merely a suite of interaction rules that must be followed once the system is stimulated externally. Thus, we must additionally define initial and external boundary conditions. Many LBM users like to invoke forcing functions. These represent body forces and may be appropriate for gravity-driven flows. This work made use of an externally controlled pressure driving force in fluid propagation studies. In the LBM, this is achieved by setting a lattice Boltzmann density gradient. 2.2. Numerical experiments 2.2.1. Spontaneous imbibition In order to probe for process sensitivity to applied wetting conditions, spontaneous imbibition experiments were performed. Without gravity present, capillaries were filled with nonwetting fluid of prescribed
Fig. 3. Comparison of theoretical and simulated permeability for straight and sinusoidally constricted capillaries: (a) permeability and aspect ratio for straight capillaries as a function of diameter, and (b) theoretical versus simulated permeability for a series of sinusoidally constricted capillaries of the same maximum diameter and wavelength, but differing minimum constriction diameter.
tion for wall density versus contact angle. The developed calibration is provided in Fig. 5 for an 11cell capillary diameter. Note the extreme sensitivity of the contact angle to the wall color density difference near F 1.8. This calibration was tested for smaller capillary sizes. The results are furnished in
Fig. 4. Examining the central longitudinal pressure profile to extract the contact angle from the LaPlace pressure jump at each fluid – fluid interface. The graphical inset shows a central cross-section of the phase map within the 3-D tube.
166
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
Fig. 5. Contact angle calibration as determined by the captive drop technique as a function of the assigned wetting density difference. Graphical inserts depict a central slice of the phase distribution within the capillary for the various points. Note that the contact angle reported is that measured through the blue (gray in these images) phase.
uniform density. At time zero, a wetting phase was allowed to contact the capillary end. Co-current imbibition of wetting phase was monitored with time.
2.2.2. Hysteresis experiments In hysteresis studies, captive drops were forced to propagate under the influence of a defined external density gradient. Wettability parameters were chosen based upon the equilibrium contact angle calibration. For propagating drops, the same pressure jump technique was used although a viscous pressure drop was superimposed. Advancing and receding angles were analyzed. Tests were conducted for different density gradients to probe for velocity effects. Similar tests were performed in sinusoidally constricted capillaries.
3. Results and discussion 3.1. Spontaneous imbibition in cylindrical capillaries
Fig. 6. Sensitivity of wettability assignment with respect to capillary size. Contact angles measured through the blue phase.
Spontaneous imbibition was observed to be sensitive to the equilibrium contact angle assigned to the walls. This sensitivity persisted beyond the point of perfect wetting. A plot of spontaneous imbibition into cylindrical capillaries is given in Fig. 7 as a function
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
167
Fig. 7. Sensitivity to wettability as indicated by spontaneous imbibition rate into a cylindrical capillary. Note that sensitivity to the wetting strength still occurs below the point of perfect wetting (h = 0). Contact angles shown are measured through the blue phase.
of wettability assignment. Although not compared directly to theoretical or experimental imbibition rates, the demonstrated sensitivity is in agreement with expectations.
faces does not reflect the same time or length scale for hydrodynamic interactions. 3.3. Propagation in sinusoidally constricted capillaries
3.2. Propagation in cylindrical capillaries Fig. 8 shows a plot of pressure profile for propagating captive drops in a cylindrical capillary for different driving forces, and thus velocities. Mass velocities were computed and used to distinguish curves in the legend. Within statistical accuracy, there was no difference observed between advancing and trailing contact angles, nor was evidence of velocity dependence displayed. These experiments were performed in a capillary of only 5 lattice units in diameter. With the walls in such close proximity, the drop shape must have been strongly influenced by conditions imposed at the walls. Still, the inability for hydrodynamics to affect the interfacial curvature is a signal that perhaps the rigor in which the equilibrium condition is applied to construct and preserve inter-
Further tests were performed for more complex, sinusoidally constricted capillaries. Without straight walls, it was anticipated to be a composite effect of equilibrium contact angle applied locally and the oscillating slope of the solid wall with respect to the tube axis. For fixed maximum diameter, wavelength, and equilibrium contact angle, propagating slugs were monitored for two different choices of minimum constriction diameter. Results are shown in Figs. 9 and 10. Clearly, interfacial motion for advancing and trailing edges becomes correlated, although the results should vary widely with slug size relative to medium wavelength. The apparent contact angles at both advancing and receding interfaces take on a range of values, challenging the meaning of apparent contact angle for convergent –divergent porous media. In both cases,
168
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
Fig. 8. Examining the central longitudinal pressure profile to extract advancing and receding contact angles, with respect to the blue phase, from captive drops propagating down a cylindrical capillary tube. No statistically significant contact angle hysteresis was observed, nor was velocity dependence evidenced over the range studied. Wetting parameter set at 2.5 to yield an equilibrium angle of 37j.
Fig. 9. Examining central longitudinal pressure jumps at leading and trailing interfaces as a function of time for a captive drop propagation down a sinusoidally constricted tube. Capillary parameters are as indicated. Wetting set for he = 65j.
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
169
Fig. 10. Examining central longitudinal pressure jumps at leading and trailing interfaces as a function of time for captive drop propagation down a sinusoidally constricted tube. Minimum capillary diameter reduced to 3 lattice units. Wetting set for he = 65j.
the maximum pressure difference is observed when the equilibrium angle is found macroscopically. Both advancing and receding angles are found to achieve this maximum. In Fig. 9, the trailing interface is even observed to take on negative pressure differences, indicating apparent wetting behavior counter to the naturally imposed wettability. This behavior could be easily predicted for capillary-dominated flow conditions as a hemispherical cap oscillates in curvature due to the combination of local equilibrium contact angle and wall curvature.
4. Conclusions 4.1. General observations The lattice Boltzmann method offers unique testing opportunities with a higher degree of control and analysis capability than physical experiments on wetting behavior. It was demonstrated that the Gunstensen method can be calibrated to yield local contact angle assignments. It also was observed that different
wettability assignments produced different kinetic behavior for numerical spontaneous imbibition experiments. 4.2. Hysteresis in LBM 4.2.1. Observations Hysteresis effects were briefly examined by simulating a propagating captive drop in a virtual capillary tube of defined equilibrium contact angle. Analysis of experiments readily yielded accurate estimates of dynamic contact angle through the use of the pressure profile across the interfaces. In the limited testing performed, no statistically significant differences were observed between advancing and receding contact angles in straight tubes. In addition, no velocity dependence of contact angles was observed. 4.2.2. Recommendations The inability to make hysteresis and velocitydependent observations challenges the current mode of interfacial orientation determination and preservation in the LBM. Although seemingly rigorous in its
170
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171
approach, the Gunstensen method performed every time step may not allow enough interplay between interfacial physics and hydrodynamics. While other methods are under investigation, introduction of a relaxation parameter may allow a balancing of the time scale for surface interactions and that for hydrodynamic interactions. Angles other than the equilibrium contact angle could additionally be imposed to force simulated behavior to conform to experimental observation, but adding ad hoc calibrations is undesirable. Rather, efforts should be made to portray the underlying physics from a mathematical perspective and allow the power of a physically meaningful fluid dynamics tool to integrate the physics over time and space. The real challenges are not in conducting flows in simple geometry, but the ones that are so complex that intuition no longer gives target behavior (Hazlett et al., 1998; Vaidya et al., 1996). Additional testing with single interfaces and other capillary sizes or shapes is in order.
t F k LBM Pc
4.3. Flow in sinusoidally constricted capillaries
The authors gratefully acknowledge the support and assistance of Drs. S.Y. Chen, D.W. Grunau, and W.E. Soll in this collaborative effort. A portion of the work was conducted using the resources of the Advanced Computer Laboratory at Los Alamos National Laboratory.
Additional simulations in sinusoidally constricted tubes indicated advancing and trailing interfaces were changing curvature resulting from interplay between a locally imposed equilibrium contact angle and the slope of the walls with respect to the central axis. If wall curvature were neglected, interpretation of the pressure jump across interfaces would be interpreted as variations in apparent contact angle. This type of analysis is not new, but the interaction between capillarity and hydrodynamics to look at hysteresis in a controlled environment with LBM was novel. The observed curvature variation does not lessen the criticism of the wettability portrayal method. It does, however, illustrate that the LBM can be used to study apparent contact angle in surface roughness or heterogeneous wettability environments. Nomenclature A interfacial tension controlling parameter d Boltzmann density, mass/volume dmax maximum capillary diameter, lattice units dmin minimum capillary diameter, lattice units db Boltzmann blue phase density, mass/volume dr Boltzmann red phase density, mass/volume D capillary diameter, lattice units
S Sw ux / k r h he ha hr
surface normal medium permeability, length2 lattice Boltzmann method capillary pressure or LaPlace pressure difference wetting strength parameter, db dr at the solid surface wetting phase saturation, volume/volume Boltzmann mass velocity, lattice units/time step medium porosity, volume/volume periodic medium wavelength, l.u. interfacial tension contact angle equilibrium contact angle advancing contact angle receding contact angle
Acknowledgements
References Cain, J.B., Francis, D.W., Venter, R.D., Neumann, A.W., 1983. Dynamic contact angles on smooth and rough surfaces. J. Colloid Interface Sci. 94 (1), 123 – 130. Collet, P., De Coninck, J., Dunlop, F., Regnard, A., 1997. Dynamics of the contact line: contact angle hysteresis. Phys. Rev. Lett. 79 (19), 3704 – 3707. de Gennes, P.G., 1985. Wetting: statics and dynamics. Rev. Mod. Phys. 57 (3/1), 827 – 863. Dussan, V.E.B., 1976. The moving contact line: the slip boundary condition. J. Fluid Mech. 77, 665 – 684. Dussan, V.E.B., Davis, S., 1974. On the motion of a fluid – fluid interface along a solid surface. J. Fluid Mech. 65, 71 – 95. Frisch, U., Hasslacher, B., Pomeau, Y., 1986. Lattice – gas automata for the Navier – Stokes equation. Phys. Rev. Lett. 56 (14), 1505 – 1508. Good, R.J., 1952. A thermodynamic derivation of Wenzel’s modification of Young’s equation for contact angles; together with a theory of hysteresis. J. Am. Chem. Soc. 74, 5041 – 5042. Grunau, D.W., 1993. Lattice Methods for Modeling of Hydrodynamics. Ph.D. Thesis, Colorado State University.
R.D. Hazlett, R.N. Vaidya / Journal of Petroleum Science and Engineering 33 (2002) 161–171 Grunau, D.W., Chen, S., Eggert, K., 1993. A lattice Boltzmann model for multiphase fluid flows. Phys. Fluids A 5 (10), 2557 – 2562. Gunstensen, A.K., Rothman, D.H., Zaleski, D.H., Zanetti, G., 1991. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43 (8), 4320 – 4327. Hazlett, R.D., 1990. Fractal applications: wettability and contact angle. J. Colloid Interface Sci. 137 (2), 527 – 533. Hazlett, R.D., 1992. On surface roughness effects in wetting phenomena. J. Adhes. Sci. Technol. 6 (6), 625 – 633. Hazlett, R.D., Chen, S.Y., Soll, W.E., 1998. Wettability and rate effects on immiscible displacement: lattice Boltzmann simulation in microtomographic images of reservoir rocks. J. Pet. Sci. Eng. 20 (3/4), 167 – 175. Hoffman, R., 1975. J. Colloid Interface Sci. 50, 228. Joanny, J.F., de Gennes, P.G., 1984. A model for contact angle hysteresis. J. Chem. Phys. 81, 552 – 562. Langaas, K., 1999. Modelling of Immiscible Two-Phase Flow in Porous Media with the Binary Fluid Lattice Boltzmann Method. Ph.D. Thesis, University of Bergen, Norway. Miksis, M.J., Davis, S.H., 1994. Slip over rough and coated surfaces. J. Fluid Mech. 273, 125 – 139.
171
Oliver, J.F., Mason, S.G., 1977. Microscopic studies on rough surfaces by scanning electron microscopy. J. Colloid Interface Sci. 60, 480 – 487. Orlandini, E., Swift, M.R., Yeomans, J.M., 1995. A lattice Boltzmann model of binary-fluid mixtures. Europhys. Lett. 32 (6), 463 – 468. Neira, M.A., Payatakes, A.C., 1979. Collocation solution of creeping Newtonian flow through sinusoidal tubes. AIChE J. 25, 725 – 730. Payatakes, A.C., Dias, M.M., 1984. Immiscible microdisplacement and ganglion dynamics in porous media. Rev. Chem. Eng. 2 (2), 85 – 174. Rothman, D.H., Keller, J.M., 1988. Immiscible cellular-automaton fluids. J. Stat. Phys. 52 (3/4), 1119 – 1127. Vaidya, R.N, Hazlett, R.D., Buckles, J.J., Coles, M.E., Grunau, D.W., Chen, S.Y., Eggert, K.G., 1996. Simulation of wettability effects on displacement processes by cellular automata methods. In: Morrow, N.R. (Ed.), Proceedings 3rd International Symposium on Evaluation of Reservoir Wettability and Its Effect on Oil Recovery. University of Wyoming, Laramie, pp. 89 – 93. Wolfram, S., 1986. Cellular automaton fluids: I. Basic theory. J. Stat. Phys. 45 (3/4), 471 – 526.