Homogeneous Volume-of-Fluid (VOF) Model for Simulating the

ARTICLE pubs.acs.org/EF

Homogeneous Volume-of-Fluid (VOF) Model for Simulating the Imbibition in Porous Media Saturated by Gas Bruno A. M. Carciofi,*,† Marc Prat,‡ and Jo~ao B. Laurindo† †

Department of Chemical and Food Engineering, Federal University of Santa Catarina, Post Office Box 476, 88040-900 Florianopolis, Santa Catarina (SC), Brazil ‡ Institut de Mecanique des Fluides de Toulouse (IMFT), Universite de Toulouse, INPT, UPS, CNRS, Avenue Camille Soula, 31400 Toulouse, France ABSTRACT: Biphasic flow in porous media is of great interest in several areas, such as petroleum and oil recovery, food impregnation, hydrology, and soil science. In this work, the penetration of an external liquid into a porous medium saturated with air, caused by the macroscopic pressure difference, was studied (impregnation). For that, a porous medium saturated with air was considered and immersed in a container with a wetting liquid. In this situation, a forced imbibition of the porous space occurs, until the equilibrium between internal and external pressure is reached, after internal gas compression. A mathematical model for this impregnation process was developed for a two-dimensional (2D) homogeneous porous medium using the volume-of-fluid (VOF) method. The medium impregnation rate calculated with the model developed in the present study was compared to the values obtained from reference analytic solutions, i.e., the equations of Lucas
Washburn and Darcy’s law. An excellent agreement was found between simulation results and the reference mathematical solutions, which indicates that the dynamics of liquid invasion can be predicted if the global pressure gradient is known. The model allows for following the interface advance for classical and nonclassical geometries, using a fixed mesh, without the need of its reconstruction at each interface displacement. This procedure is timesaving and inputs robustness to the computational model.

’ INTRODUCTION Biphasic flow in porous media is important in many industrial and geological applications, such as petroleum recovery, vacuum impregnation of foods, hydrology, and soil science. In nonmiscible biphasic systems, different flow patterns are observed, depending upon the fluid wettability, the displacement velocity, and also the ratios of the viscosities and densities of the fluids. Displacements into porous media are a consequence of the pressure difference, originating from a macroscopic gradient imposed on the system or natural phenomena (capillarity and gravitational force). Nowadays, there is an increasing interest in precisely evaluating the magnitude and role of the capillary effect in infiltration processing (spontaneous capillary imbibitions).1 Two methods are commonly used to follow the interface during a biphasic flow: (i) Lagrangean-based methods that follow the interface during the displacement and (ii) Eulerian-based methods, used to describe the fluid movement from an inertial reference frame. Lagrangean methods allow for accessing displacement hydrodynamic details and following interfaces with precision but are difficult to use for situations with two-fluid snap-off and coalescence and when the interface length changes continuously.2 Eulerian methods use scalar function to determine interface position that allows for using a fixed mesh, on which the interface displacements take place.3 The socalled volume-of-fluid (VOF) method is the Eulerian method used by the computational code JADIM, developed at the Institut de Mecanique de Fluides de Toulouse (IMFT), Toulouse, France. This is a numerical approach suitable for tracking interface deformation along with direct simulation of r 2011 American Chemical Society

flow of liquid and gas phases. The code uses a coupled levelset/VOF interface capturing method to determine the volumetric fraction of one phase (liquid phase for example) by solving the transport equations for this phase, without the need for interface reconstruction. Interface reconstruction for three-dimensional problems is a complex problem.3 Classically, immiscible fluid displacements in porous media may take one of the three basic forms: viscous fingering, capillary fingering, or stable displacement. The first displacement type occurs during the fast drainage of a liquid by a gas, when the pressure drop is the dominant phenomena. Two-phase displacement patterns called capillary fingering are present at a low capillary number, when the viscous forces are negligible in both fluids and the capillarity effects dominate the process. In this case, the channel radius distribution controls the shape of the invasion front. A stable displacement pattern occurs when the much higher viscosity of the injected fluid stabilizes the invading front, leading to a piston-type displacement. At this condition, there is neither viscous fingering nor capillary fingering.4 In this work, the penetration of an external liquid into a porous medium saturated with air, caused by a coupled macroscopic pressure difference and capillary pressure, was studied. For that, a porous medium saturated with air was considered immersed in a container with a wetting liquid. In this situation, a forced imbibition of the porous solid occurs, until the equilibrium between the internal and external pressure is reached, because Received: February 14, 2011 Revised: April 5, 2011 Published: April 08, 2011 2267

dx.doi.org/10.1021/ef200233j | Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

ARTICLE

of the internal gas compression. A mathematical model for describing this impregnation process was developed for a twodimensional (2D) porous media using the VOF model, valid for two-phase displacements with high capillary number and/or high viscosity of the invading fluid. This is the case of porous media initially saturated with a gas and invaded by a liquid, leading to stable displacement patterns. The impregnation rates calculated with the model developed in the present study were compared to the values obtained from reference analytic solutions, i.e., the equations of Lucas
Washburn5,6 and Darcy’s law in radial coordinates.7

’ MATERIALS AND METHODS Mathematical Modeling. The system under analysis is idealized as an isotropic and homogeneous porous medium, with constant permeability (K) and porosity (ε) and with pore sizes represented by an average radius (Rc). It was considered that the system is isothermal and that (i) the liquid
gas interfacial tension is constant, (ii) the two fluids are Newtonian and immiscible, and their viscosities are constant, (iii) the liquid is incompressible, (iv) the gas mixture behaves ideally, and (v) the solid is completely inert. The gas displacement inside the porous media because of the liquid invasion can be described as a quasistationary state, which implies that the gas density is spatially homogeneous. This is because the pressure relaxation in the gas phase is very fast compared to the interface displacement. The physical and mathematical representation of the displacement of a gas by a liquid inside a porous media may be given by eqs 1
7. Equations 1 and 2 are the continuity equation and the momentum conservation equation,8 respectively. Equation 3 represents the equation of state to ideal gas behavior, while eq 4 is used to calculate the capillary pressure through the interface (Laplace equation9). Equation 5 represents the sum of capillary and liquid-phase pressures, called corrected pressure (PL*). The pressure gradient through the gas
liquid interface and the interface displacement velocity are eqs 6 and 7, respectively rUk ¼ 0 μ rPk þ Fk g þ k Uk ¼ 0 K PG ¼ FG

RT M

2σLG cos θ Rc

ð1Þ

At the system boundaries, the pressure source condition was considered (Dirichlet boundary condition), represented by eq 8 PL  ¼ Po þ Pc

ð8Þ

where Po is the external pressure imposed on the porous medium system. Formulation of the VOF Model. The Eulerian approach does not allow for the numeric solution of local balance equations (eqs 1
8), because the two fluids can be present in a VOF method used to perform the calculation.2 An alternative is using a homogeneous model that describes the displacement of two immiscible fluids as a single equivalent fluid. This way, the specific gravity and viscosity of the fluid, defined by the homogeneous model, vary through the interface accordingly with the physical properties of each fluid. The model used to simulate the impregnation of a porous medium by a liquid is based on the work by Benkenida,10 using a phase-indicating function and the conservation of mass and momentum at the liquid and gas phases. This model allows for evaluating the volumetric fraction of a k phase (Rk), the interface position, and the fraction of the porous medium occupied by liquid (XL). In a liquid
gas displacement, at every instant t, at every position x in the system, the phases can be identified by a phase-indicating function jk(t,x), such that jk = 1 if the k phase is present at the position x, at the instant t. Otherwise, jk = 0. Because jk is relative to the fluid, its derivative following the interface displacement is null,10 leading to Djk þ UΓ rjk ¼ 0 Dt From eqs 1, 7, and 9 ε

Djk þ rðjk Uk Þ ¼ 0 Dt

ð3Þ

ð10Þ

For biphasic displacements, one can write DjG Dj ¼
L Dt Dt

ð11Þ

Hence, the summation of eq 10 applied to both phases leads to rðjG UG þ jL UL Þ ¼ 0

ð12Þ

The summation of eq 9 applied to both phases and using eq 11 leads to rjG ¼
rjL

ð2Þ

ð9Þ

ð13Þ

From the momentum equation (eq 2) applied to both phases, multiplied by the respective values of jk, and from eqs 6 and 13, one can write ! 1 jk Pk þ g Fk jk þ μk Uk jk ¼ 0 ð14Þ r K k ¼ L, G k ¼ L, G k ¼ L, G

∑

∑

∑

ð4Þ

PL  ¼ Pc þ PL

To write the homogeneous model, a filter must be applied to the previous equations. The filter has a Y operator applied around the x position, such that10

ð5Þ

Y ðx
x 0 Þ g 0, " x and " x 0

ðPG
PL ÞΓ ¼ 0

ð6Þ

Pc ¼

UG UL n¼ n ¼ UΓ n ε ε

ð7Þ

where U, P, F, and μ are the velocity vector, pressure, specific gravity, and dynamic viscosity, respectively, defined for liquid (L) and gas (G) phases or at the gas
liquid interface (Γ), g is the gravity acceleration, R is the universal gas constant, T is the system temperature, M is the average molar mass of the gas phase, σLG is the interfacial tension at the liquid
gas interface, θ is the contact angle, and n is the normal vector to the liquid interface.

Therefore

Z ϑ

Yðx
x 0 Þdx 0 ¼ 1

ð15Þ

ð16Þ

where ϑ is the filter characteristic volume, which is a dimension of a calculus element used for solving equations numerically. The filtered expression of the function jk is given by Z Yðx
x 0 Þjk ðt, x 0 Þdx 0 ¼ Æjk æðt, xÞ ¼ Rk ðt, xÞ ð17Þ ϑ

The new function Rk(t,x), known as the presence fraction of the k phase, is the instantaneous volumetric fraction of the k phase, on a radius x0 2268

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

ARTICLE

Figure 1. Comparison of numerical and Lucas
Washburn solutions (one-dimensional problem).

Table 1. Parameters Used in the Simulations Related to Figures 1 and 2 parameter

value

K

1.25  10
11

unit m2

PL*

1.01325  10

PGo

3.3  103

Pa

μG

1.8  10
5

Pa s

μL

1.07

Pa s

FL Rc

1  103 1  10
5

kg m
3 m

ε

2  10
1

fraction

dimension of the porous media

9.95  10
3

m

5

Pa

Figure 3. RG evolution for a 2D displacement in a disk (symmetric problem). (a) t = 0; (b) 10 iterations, t = 2.7  10
4 s; (c) 50 iterations, t = 2.0  10
2 s; (d) 500 iterations, t = 2.1  10
1 s; (e) 5000 iterations, t = 2.1 s; (f) 10 000 iterations, t = 4.2 s; (g) 13 000 iterations, t = 5.4 s; and (h) 19 500 iterations, t = 8.0 s.

Figure 2. Influence of numerical refining on the solution of the onedimensional problem, using grids with 100 and 200 equally spaced divisions. Comparison to the solution of Lucas
Washburn.

ϑ

Y ðx
x 0 Þjk ðt, x 0 ÞΦk ðt, x 0 Þdx 0 ¼ Rk ðt, x 0 ÞÆΦk æðt, xÞ

ð19Þ

ÆPæ ¼ RG ÆPG æ þ RL ÆPL æ

ð20Þ

Æ μæ ¼ RG μG þ RL μL

ð21Þ

ÆFæ ¼ RG ÆFG æ þ RL FL

ð22Þ

Applying this filter to the equations describing the problem leads to a new set of equations rÆUæ ¼ 0

around the position x. For a hypothetical variable Φk(t,x0 ), eq 18 represents the filtered value of this variable, ÆΦkæ(t,x), called the local averaged value. Z

ÆUæ ¼ RG ÆUG æ þ RL ÆUL æ

rÆPæ þ ÆFæg þ

ε

ð18Þ

The velocity, pressure, viscosity, and specific gravity used in the homogeneous model (local averaged values) are given by eqs 19
22.

1 ÆμæÆUæ ¼ 0 K

DRG þ ÆUærRG ¼ 0 Dt ÆFG æ ¼ ÆPG æ

2269

_ M RT

ð23Þ ð24Þ

ð25Þ

ð26Þ

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

ARTICLE

Table 2. Parameters Used in the Simulations Related to Figures 3 and 4 parameter

value

unit
11

K

1.25  10

PL* PGo

1.01325  10 3.3  103

Pa Pa

μG

1.8  10
5

Pa s

μL

1.07

Pa s

FL

1  103

kg m
3

Rc

1  10
5

m

ε

2  10
1

fraction

ro

9.83  10
3

m

m2 5

Figure 5. Evolution of XL for a 2D displacement in a disk (symmetric problem). Comparison between numerical solution (400  400 grid) and the reference solution given by eq 30. Figure 4. Pressure (in Pa) evolution for a 2D displacement in a disk (symmetric problem). (a) t = 0; (b) 10 iterations, t = 2.7  10
4 s; (c) 50 iterations, t = 2.0  10
2 s; (d) 500 iterations, t = 2.1  10
1 s; (e) 5000 iterations, t = 2.1 s; (f) 10 000 iterations, t = 4.2 s; (g) 13 000 iterations, t = 5.4 s; and (h) 19 500 iterations, t = 8.0 s. Z ÆPG æ ¼ ÆPGo æ ZV

RGo dV RG dV

ð27Þ

V

where PGo and RGo are the pressure and gas volumetric fraction at the initial instant, while V is the global volume of the porous media (system volume). The volumetric fraction of liquid that impregnates the porous media is given by eq 28. 0 Z 1 R dV G B C V C ð28Þ XL ¼ εB @1
V A

Reference Equations for Imbibition of Porous Media. Reference solutions can be used for validating the homogeneous model. For the 1D case, the equation of Lucas
Washburn5,6 was preferred, while for the 2D case, an equation derived from the Darcy equation was chosen.7 Lucas
Washburn equation (1D case): The equation proposed by Lucas5 and Washburn6 is derived from a force balance between the

Figure 6. Evolution of XL for a 2D displacement in a disk (symmetric problem). Influence of grid refining (60  60, 200  200, and 400  400) on the numerical solution and comparison to the reference solution. forces that drive the fluid displacement (capillary forces) and the resisting forces (viscous forces), under laminar flow. Gravity forces can drive or hinder fluid displacement, depending upon the direction of the displacement. The displacement of a gas in a cylindrical capillary by a liquid (with much higher viscosity than gas viscosity) is given by eq 29 8μL z 2270

dz ¼ 2Rc σLG cos θ ( Fe dt

ð29Þ

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

ARTICLE

Figure 7. Evolution of RG for a 2D non-classical geometry at different steps: t = 0, 0.24, 1.0, 2.0, 5.0, 7.0, 10.0, 13.0, 17.0, 20.0, and 30.0 s, respectively. where Rc is the capillary radius, z is the instantaneous position of the interface inside the capillary, and Fe is an external force acting on the system (e.g., gravitational force). Radial imbibition in a porous disk (2D case): Darcy’s law was used to describe the radial (2D) imbibition of a porous disk initially saturated with air, leading to an equation that describes the instantaneous position of the gas
liquid interface (eq 30) "  #
1 dr K r ¼ ðPL, ro
PL, r Þ ln ð30Þ dt μL rε ro where r is the radius of the interface position (i.e., the position of the imbibition front), ro is the radius of the interface at t = 0, and PL,r and PL,ro are the liquid pressures relative to r and ro. Numerical Solution. Numerical solutions applied to the homogeneous model should satisfy at least these four conditions: 10 (i) stability, (ii) conservation of RG that leads to local and global mass conservation (incompressible liquid and pseudo-stationary condition for the gas phase), (iii) positivity of the RG function, and (iv) spatial precision, which reduces the numerical diffusion caused by discontinuity of the RG function. The solution of the homogeneous VOF model was performed using the JADIM code, developed at the IMFT, Toulouse, France. The flux-corrected transport (FCT) methodology is used by the JADIM code. Initially proposed by Boris and Book,11 the FCT numerical scheme was developed to solve collision problems and was modified by Zalesak12 for treating numerical diffusion. Benkenida11 and Bonometti2 reported works about the improvement of the solution of the phase presence fraction using the FCT scheme. In the JADIM code, the numerical solution for RG was obtained using a staggered grid. The discretization of the presence fraction equation (eq 25) was performed by a second-order finite-volume method.3,13 The pressure field was calculated numerically by discretizing eqs 23 and 24 by an implicit finitevolume method.14 Statistical Parameters. The statistical parameters mean square error (MSE) (eq 31) and bias factor (eq 32) were used for comparing

values calculated from the reference solutions and from the homogeneous model and for comparing different numerical solutions.

∑ðreference value
simulated valueÞ2

ð31Þ

bias ¼ 10∑logðreference value=simulated valueÞ=number of data

ð32Þ

MSE ¼

number of data

’ RESULTS AND DISCUSSION Simulation results of impregnation of a one-dimensional porous medium, initially saturated with gas and open at only one side, in contact with the impregnating liquid at constant pressure, are given in Figure 1. The value of the fraction, XL, calculated by the homogeneous VOF model, using a computational grid of 100 equally spaced divisions, is compared to the value given by the equation of Lucas
Washburn. The parameters used in the simulations are given in Table 1, representing a vacuum impregnation. Figure 1 shows that the homogeneous model results agree very well with the solution of Lucas
Washburn. A grid sensitivity analysis is shown in Figure 2, comparing the reference solution to numerical solutions using grids with 100 and 200 equally spaced divisions. The grid refining improved the numerical solution, leading to a 2.6-fold increase in the computation time in an Intel Xeon 3.06 GHz personal computer (PC), on Linux operational system (OS). For simulations with 60 000 iterations, the MSE value varied from 1.81  10
6 to 1.25  10
6, while the bias factor did not vary significantly (from 1.008 to 1.009) when the grid was changed from 100 to 200 divisions. The simulation results on the impregnation dynamics of a two-dimensional porous medium (porous disk, represented by 1 /4 of the disk, because of the geometric symmetry) is given by the evolution of RG during the impregnation process and is shown in Figure 3, while the pressure fields are shown in Figure 4. 2271

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

ARTICLE

Figure 8. Pressure (in Pa) evolution for a 2D non-classical geometry at different steps: t = 0, 0.24, 1.0, 2.0, 5.0, 7.0, 10.0, 13.0, 17.0, 20.0, and 30.0 s, respectively.

Table 3. Parameters Used in the Simulations Related to Figures 7 and 8 parameter

value

unit

K

1.25  10
11

m2

PL*

1.01325  105

Pa

PGo μG

3.3  103 1.8  10
5

Pa Pa s

μL

1.07

Pa s

FL

1  103

kg m
3

Rc

1  10
5

m

ε

2  10
1

fraction

In this case, the system had four open sides accessible to the impregnating liquid. The simulation parameters used in these simulations are given in Table 2. All simulations were performed using a 400  400 grid. Figure 3 shows that the liquid
gas interface is stable and well-defined, indicating that the model and numeric solution used to solve the problem are appropriate for this kind of phenomenon. The pressure field shown in Figure 4 is coherent, in agreement with the displacement geometry. The results of XL achieved with the homogeneous VOF (400  400 grid) is in agreement with the reference solution given by eq 30, as shown in Figure 5. A grid sensitivity analysis is shown in Figure 6, comparing the reference solution to numerical solutions using 60  60, 200  200, and 400  400 grids. The numerical results are in agreement with the solution represented by eq 30, and the agreement is better when the grid is refined. Figures 7 and 8 show the evolutions of RG and of the pressure field during the impregnation of a 2D porous medium with a nonclassical geometry, for which there is no analytical reference solution. As in previous simulation, the system was fully open and the liquid could enter by all sides. The parameters used in these

simulations are given in Table 3. Figure 7 shows that the homogeneous/VOF model allowed for following the interface displacement during time, which was stable and well-defined during the simulation. The results show that the geometry of the region occupied by the gas phase at the end of the simulation is different from the geometry occupied by the gas phase at the beginning of the simulation (initial condition). The pressure field shown in Figure 8 is coherent, in agreement with the displacement geometry.

’ CONCLUSION The agreement between the results obtained with the homogeneous VOF model and the values found with reference mathematical solutions validates its use for simulating the impregnation of one-dimensional or two-dimensional porous media by a liquid, driven by global pressure gradients. It is the case of vacuum impregnation of porous media. The model also allows for simulating vacuum impregnation in non-classical geometries using a fixed grid; i.e., it is not necessary to reconstruct the grid at each interface displacement, leading to computational time saving and code robustness. The homogeneous model allows for following very well the gas
liquid interface during the impregnation, which is really important for many situations of impregnations of interest. We have only considered two-dimensional examples. It is obvious that the method can be straightforwardly extended to compute the impregnation process in porous medium of complex three-dimensional shapes. ’ AUTHOR INFORMATION Corresponding Author

*Fax: 55-48-3721-9687. E-mail: [email protected]. 2272

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273

Energy & Fuels

’ ACKNOWLEDGMENT The authors thank Dominique Legendre and Annaïg Pedrono (IMFT, Toulouse, France) for their useful help on the use of the JADIM code and CAPES, CNPq, Brazil for financial support. ’ NOMENCLATURE 1D = one-dimensional 2D = two-dimensional FCT = flux-corrected transport G = gas phase L = liquid phase MSE = mean square error OS = operational system VOF = volume-of-fluid Γ = interface Æ æ = average value Fe = external force acting on the system (N) g = gravity acceleration (m s
2) K = permeability (m2) M = average molar mass of the gas phase (kg mol
1) n = normal vector to the liquid interface (nondimensional) P = pressure (Pa) Pc = capillary pressure (Pa) PG = pressure at the gas phase (Pa) PGo = pressure at initial time (Pa) Pk = pressure at the k phase, where k = L or G (L, liquid; G, gas) (Pa) PL = pressure at the liquid phase (Pa) PL* = corrected pressure at the liquid phase (Pa) Po = external pressure (Pa) r = radius of the interface position (m) ro = radius of the interface at initial time (m) R = universal gas constant (J mol
1 K
1) Rc = average pore size radius (m) t = time (s) T = system temperature (K) U = velocity vector (m s
1) UG = velocity vector for the gas phase (m s
1) Uk = velocity vector for the k phase, where k = L or G (L, liquid; G, gas) (m s
1) UL = velocity vector for the liquid phase (m s
1) UΓ = velocity vector for the interface (m s
1) V = global volume of the porous medium (m3) x = position vector (m) XL = fraction of the porous medium occupied by liquid (fraction) Y = filter operator (nondimensional) z = instantaneous position of the interface inside the capillary (m) RG = volumetric fraction of the gas phase (fraction) RGo = volumetric fraction of the gas phase at initial time (fraction) Rk = volumetric fraction of the k phase, where k = L or G (L, liquid; G, gas) (fraction) RL = volumetric fraction of the liquid phase (fraction) ε = porosity (fraction) θ = contact angle (rad) μ = dynamic viscosity (Pa s) μG = dynamic viscosity for the gas phase (Pa s) μk = dynamic viscosity for the k phase, where k = L or G (L, liquid; G, gas) (Pa s) μL = dynamic viscosity for the liquid phase (Pa s)

ARTICLE

F = specific gravity (kg m
3) FG = specific gravity for the gas phase (kg m
3) Fk = specific gravity for the k phase, where k = L or G (L, liquid; G, gas) (kg m
3) FL = specific gravity for the liquid phase (kg m
3) σLG = interfacial tension at the liquid
gas interface (N m
1) jG = gas-phase-indicating function (nondimensional) jk = k-phase-indicating function, where k = L or G (L, liquid; G, gas) (nondimensional) jL = liquid-phase-indicating function (nondimensional) Φk = hypothetical variable, where k = L or G (L, liquid; G, gas) (nondimensional) ϑ = filter characteristic volume (m3)

’ REFERENCES (1) Cai, J.; Yu, B.; Zou, M.; Luo, L. Energy Fuels 2010, 24, 1860–1867. (2) Bonometti, T. Developpement d’une methode de simulation d’ecoulements a bulles et a gouttes. Ph.D. Thesis, Institut National Polytechnique de Toulouse, Toulouse, France, 2005. (3) Bonometti, T.; Magnaudet, J. Int. J. Multiphase Flow 2007, 33 109–133. (4) Lenormand, R.; Touboul, E.; Zarcone, C. J. Fluid Mech. 1988, 189, 165–187. (5) Lucas, R. Kolloid-Z. 1918, 23, 15–22. (6) Washburn, E. W. Phys. Rev. 1921, 17, 273–283. (7) Hyv€aluoma, J.; Raiskinm€aki, P.; J€asberg, A.; Koponen, A.; Kataja, M.; Timonen, J. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2006, 73, No. 036705. (8) Sisson, L. E.; Pitts, D. R. Elements of Transport Phenomena; McGraw-Hill: New York, 1972. (9) De Gennes, P.-G.; Brochard-Wyart, F.; Quere, D. Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves; Springer-Verlag: New York, 2003. (10) Benkenida, A. Developpement et validation d’une methode de simulation d’ecoulements diphasiques sans reconstruction d’interface— Application a la dynamique des bulles de Taylor. Ph.D. Thesis, Institut National Polytechnique de Toulouse, Toulouse, France, 1999. (11) Boris, J. P.; Book, D. L. J. Comput. Phys. 1973, 18, 38–69. (12) Zalesak, S. T. J. Comput. Phys. 1979, 31, 335–362. (13) Magnaudet, J.; Rivero, M.; Fabre, J. J. Fluid Mech. 1995, 284, 97–135. (14) Patankar, S. V. Numerical Heat Transfer and Fluid Flow; HPC (Taylor and Francis Group): Saint Leonards-on-Sea, U.K., 1980.

2273

dx.doi.org/10.1021/ef200233j |Energy Fuels 2011, 25, 2267–2273