Ind. Eng. Chem. Res. 2006, 45, 3311-3319
3311
Transient Natural Convection Induced by Gas Diffusion in Liquid-Saturated Vertical Porous Columns Zhaowen Li, Mingzhe Dong,* and Ezeddin Shirif Faculty of Engineering, UniVersity of Regina, Regina, Saskatchewan, S4S 0A2, Canada
This paper focuses on the development of natural convection mass transfer induced by gas diffusion in a liquid-saturated vertical porous column. The porous column is confined with two impermeable end faces and an isosolutal side surface. Mass transfer by natural convection is studied by numerically solving the conservation equations of mass, momentum, and diffusing species. To describe the transient behavior of natural convection, transient and integrated Sherwood numbers, i.e., Sht and Shi, are defined on the basis of a parallel pure diffusion case. The evolution of the natural convection process is systematically studied for the aspect ratios ranging from 0.125 to 8 and Rayleigh numbers from 10 to 1000. It is found that the Sherwood number increases with increasing Rayleigh number, whereas it increases first and then decreases with increasing aspect ratio. The results confirm the assumption in a previous study that the effect of natural convection is negligible on the diffusion process in the gas effective diffusion coefficient measurement,1 and also provide important implications on how to prepare core samples in laboratory for gas effective diffusion coefficient measurements. Introduction The dissolution of a gas in a liquid usually causes a density change of the liquid phase, which may result in natural convection when a gas diffuses into a liquid-saturated porous medium. The natural convection mass transfer of gases in liquidsaturated porous media has many important applications in various processes of petroleum, environmental, and chemical engineering. The phenomenon of natural convection, primarily natural convection heat transfer, in liquid-saturated porous media has been widely investigated.2-11 The theories and results for natural convective heat transfer can be well-adapted to analyze natural convection mass transfer due to the similarities they bear;12 however, the behavior of transient natural convection heat transfer in a vertical, cylindrical, liquid-saturated porous medium remains uninvestigated, even though the steady-state heat convection for such configurations has been studied.8,9 Therefore, from the point of view of theoretical research, the study of transient natural convection for either heat or mass transfer in vertical cylindrical liquid-saturated porous media is needed to understand this fundamental phenomenon. Moreover, one immediate application of such a study is found in the measurement of the gas effective diffusion coefficient in liquidsaturated porous media,1 in which the effect of natural convection on the measured results needs to be examined. The effective diffusion coefficient of gas in liquid-saturated porous media (Deff) under high pressures is a requisite parameter for predicting the gas migrations and distributions in geological formations in oil and gas engineering or carbon dioxide (CO2) geological sequestration practices. Despite that, the method has not been well-developed for the measurement of Deff under highpressure conditions. To determine Deff under high pressures, a method was presented in a recent study,1 in which a liquidsaturated porous rock column with two sealed end faces was used as the test sample so that the gas diffused into the rock column only through the side surface. In the measurement, the * To whom correspondence should be addressed. Tel.: +1(306)337-2269. Fax: +1(306)585-4855. E-mail:
[email protected].
column was housed in a diffusion cell filled the high-pressure test gas, CO2, and Deff was determined from the recorded gas pressure decay history. That method is of great potential for further applications in (1) measuring Deff under high pressures for different gas/liquid systems; (2) determining the diffusive tortuosity factors of the porous rock by comparing the diffusion coefficient, D, of a gas in a liquid and the Deff in the liquidsaturated rock sample; and (3) determining CO2 diffusion coefficient in bulk liquid phase at high pressures, since CO2 diffusion coefficients in most bulk liquids are difficult to measure by the conventional method due to the interference of natural convection caused by density increase as CO2 dissolves.13-15 Therefore, for further study of gas diffusion in liquid-saturated porous media or CO2 diffusion in liquid phase through the previously presented method,1 a quantitative examination of natural convection on the measured diffusion coefficients and solutions to minimize the convection effect are of practical importance. The objective of the present study is to numerically investigate the transient natural convection in a liquid-saturated vertical porous column induced by the solute diffusing into the column through the side surface area. The effect of natural convection on the diffusion process is quantitatively examined. Theoretically, the evolution of the natural convection process induced by solute diffusion in a vertical porous column can be thoroughly understood. Mathematical Formulation The problem considered herein is the natural convection mass transfer induced by the diffusion of a gas into a liquid-saturated vertical porous column. The configuration and the coordinates of the porous column are given in Figure 1, where the two end faces are sealed and impermeable to both diffusion and fluid flow. When the solute concentration at the side surface is set at a constant level, the diffusion of the solute into the liquidsaturated porous column occurs, and consequently, natural convection is generated if the solution density changes as the solute dissolves in the liquid saturating the porous medium. The equations governing the conservation of mass, momentum, and
10.1021/ie0510829 CCC: $33.50 © 2006 American Chemical Society Published on Web 03/25/2006
3312
Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006
yield a momentum equation that determines the coupling of the concentration field and the flow field, written in scalar form as
∂u′ ∂V′ F0Kgβ ∂c′ ) . ∂z′ ∂r′ µ ∂r′
(5)
Substitution of eq 1 into eq 3 yields, in scalar form,
∂c′ ∂c′ ∂c′ + u′ + V′ ) Deff∇2c′ ∂t ∂r′ ∂z′
(6)
After the above manipulations, the governing equations, eqs 1-3, are reduced in scalar form to eqs 5 and 6. For the cylindrical system shown in Figure 1, the boundary conditions are given as
u′ ) 0, c′ ) c′0, at r′ ) r′0 u′ ) 0, V′ ) 0, Figure 1. Schematic of liquid-saturated porous column subjected to a constant solute concentration at side surface.
diffusing species in the above process are given in vector form as
∂F + ∇‚(Fv′) ) 0 ∂t
(1)
K v′ ) - (∇P + Fg) µ
(2)
∂c′ ) Deff∇2c′ - ∇‚(c′v′) ∂t
(3)
∂c′ ) 0, at z′ ) 0 and z′ ) H′ ∂z′
(7)
where c′0 is the concentration at the surface of the porous column. Initially, the liquid contained in the porous medium is motionless, and the solute saturation in the liquid saturating the porous medium is c′i. At t ) 0, the concentration at the side surface of the column is suddenly raised to a constant value c′0. Thus, the initial conditions are written as
u′ ) V′ ) 0, c′ ) c′i for r′ < r′0
where, F is the liquid-phase density, t is the time, v′(u′,V′) is the velocity vector (u′ and V′ are the horizontal and vertical velocity components), K is the permeability of the porous medium, µ is the liquid viscosity, P is the pressure, g is the gravitational acceleration, c′ is the concentration, and Deff is the effective diffusion coefficient of the solute in the liquidsaturated porous medium. In the present study, the effective diffusion coefficient, Deff, is assumed to be a constant. This assumption applies for the dilute systems under constant pressure and temperature conditions.12 The constituent equation takes a linear form as
F ) F0[1 + β(c′ - c′i)]
∂c′ ) 0, at r′ ) 0 ∂r′
u′ ) V′ ) 0, c′ ) c′0 for r′ ) r′0 at t ) 0 Numerical Solution
To solve the problem described in the preceding section, the following dimensionless variables are defined
τ)
r′ z′ u′ V′ t , r) , z) , u) , V) , r′ r′ D /r′ D r′0 /Deff 0 0 eff 0 eff/r′0 c′ - c′i (9) c) c′0 - c′i 2
where τ, r, z, u, and V are the dimensionless variables corresponding to t, r′, z′, u′, and V′, respectively. Introducing the stream function, ψ,
u)-
(4)
where c′i is the reference concentration (the initial solute concentration in the liquid contained in the porous medium), F0 is the density corresponding to the concentration, c′i, and β is the volumetric concentration expansion coefficient. β > 0 and β < 0 indicate a density increase and decrease, respectively, due to the solute dissolution, with the former being the interest of this study. Invoking the Boussinessq approximation16 that the density variations with concentration in all the terms are neglected except in those that generate fluid motions (the buoyancy term Fg in eq 2), eqs 1-4 can be simplified though the following operations. Similar to the procedure in the literature,3,4 eliminating the pressure term P in eq 2 and then combining with eq 4
(8)
1 ∂ψ 1 ∂ψ , V) r ∂z r ∂r
(10)
the governing equations, eqs 5 and 6, have the following dimensionless form, respectively,
(∇ - 2r ∂r∂ )ψ ) Ra r∂c∂r
(11)
∂c ∂c ∂c + u + V ) ∇2c ∂τ ∂r ∂z
(12)
2
and
where Ra is the Darcy-modified Rayleigh number based on the radius of the porous column,
Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3313
Ra )
gKF0β(c′0 - c′i)r′0 µDeff
(13)
Table 1. Grid Fineness and Time Steps Used in Numerical Calculations H 4 to 8
The boundary conditions, given by eq 7, are reduced correspondingly to
1/8 to 2
c ) 1 at r ) 1 ∂c ) 0 at r ) 0 ∂r ∂c ) 0 at z ) 0 and z ) H ∂z ψ ) 0 at all the boundaries, r ) 1, z ) 0, and z ) H (14) where H ) H′/r′0 is the aspect ratio of the system. The initial conditions given in eq 8 are reduced to
ψ ) 0, c ) 0 for r < 1 ψ ) 0, c ) 1 for r ) 1 at τ ) 0
(15)
The formulation of the dimensionless form of the problem (eqs 11-15) shows that the solution of eqs 11 and 12 depends on two dimensionless variables, Ra and H, which are systematically examined in the next section. Similar to the practices in the literature,3,6,8,16 the governing partial differential equations, eqs 11 and 12, were solved by a finite difference method (FDM). A two-step alternating direction implicit (ADI) method was used in solving the species conservation equation, eq 12, and the point successive over relaxation method (PSOR) method was used to calculate the stream function ψ from eq 11. Central difference approximations, which have second-order accuracy, were used for the internal grid points in the discretization of eqs 11 and 12, whereas the backward difference was used at the boundaries r ) 0 and r ) 1. Forward difference was used in discretizing the time τ. For each time step, the conservation equation, eq 12, was first solved using the previous time step velocities to obtain the concentration distributions, and then, the stream function ψ was calculated from the momentum equation, eq 11, with the updated concentrations. Once ψ was obtained, the secondary velocities u and V were updated with the new values of ψ by eq 10 and, subsequently, used in updating the concentration field by eq 12. The iterative numerical calculation for each time step was repeated until the following convergence criteria were satisfied, k+1
ψi,j
|
- ψi,j
k
ψi,j
k+1
ci,j |max e and |
- ci,j
k+1
k
ci,j
k+1
|max e (16)
where i and j are the indexes of the mesh points in r and z direction, respectively, and k stands for the iteration cycles. is the prescribed tolerance of the relative errors between two iteration cycles, which was set as 10-4 in this study. Owing to the symmetry of the problem, only half of the twodimensional domain, i.e., 0 e r e 1, 0 e z e H, was of interest in the development of numerical solution. This area was further divided into I × J grids, where I and J are the maximum indexes of the mesh points in r and z direction, respectively. The fineness of the grid depended on the Rayleigh number, Ra, and the aspect ratio H. As the Rayleigh number increases, a finer grid system and smaller time steps were used to provide sufficient accuracy. The grid configuration and the time steps used in all the calculations of this study are summarized in Table 1.
Ra
I×J
∆τ
10 to 100 200 to 600 1000 10 to 100 200 to 600 1000
61 × 101 81 × 161 101 × 201 61 × 51 81 × 81 101 × 101
0.000 72 0.000 48 0.000 24 0.000 72 0.000 48 0.000 24
Analogous to Nusselt numbers in heat transfer, the Sherwood number is used to describe the intensity of mass transfer due to natural convection. For most of the natural convection heat transfer studies,10 constant temperature differences crossing the porous layers, enclosures, or annuli persist during the natural convection process such that the Nusselt numbers can be defined on the basis of the pure conduction under the corresponding temperature differences. In this study, however, the mass transfer into a liquid-saturated porous column through both pure diffusion and natural convection are transient processes. As a result, a constant concentration difference, which can be used as a reference to define the Sherwood number for the system, does not exist. Thus, in this study, the Sherwood number is defined on the basis of a parallel pure diffusion process (an assumed base case in which the pure diffusion starts at the same time as that of the natural convection case). The transient Sherwood number Sht and the integrated Sherwood number Shi are defined, respectively, to describe the transient and overall intensity of a natural convection process. The transient Sherwood number is given as
[∫ ( ) ] [∫ ( ) ] ∂c dz ∂r r)r0 H ∂c dz 0 ∂r r)r 0 H
nc Sht ) ) nd
0
c
(17)
d
where nc and nd are the mass transfer rates into the porous column through the side surface area at time τ through natural convection and the parallel pure diffusion process, respectively. The integrated Sherwood number Shi is defined in derivative or integral form, respectively, as
[∫ ∫ (∂c∂r) [∫ ∫ (∂c∂r)
], dz dt]
(18a)
∫0H ∫01 2πrc dr dz]c ∫0H ∫01 2πrc dr dz]d
(18b)
τ
qc Shi ) ) qd
0
τ
0
H
0
r)1
H
0
r)1
dz dt
c
d
and
qc [ Shi ) ) qd [
where qc and qd are the accumulated masses of the solute entered into the porous column at time τ through natural convection and the parallel pure diffusion process, respectively. For the systems free of convection (only diffusion occurs), both Sht and Shi have the value of unity. In this study, the values of Shi obtained from eqs 18a and 18b through the numerical computation agreed well with each other, with relative differences of
