Transp Porous Med (2013) 97:395–408 DOI 10.1007/s11242-013-0130-5

Onset of Buoyancy-Driven Convection in a Fluid-Saturated Porous Medium Bounded by a Long Cylinder Min Chan Kim

Received: 8 October 2012 / Accepted: 21 January 2013 / Published online: 7 February 2013 © Springer Science+Business Media Dordrecht 2013

Abstract A theoretical analysis of convective instability driven by buoyancy forces under the transient concentration fields is conducted in an initially quiescent, liquid-saturated, cylindrical porous column. Darcy’s law and Boussinesq approximation are used to explain the characteristics of fluid motion and linear stability theory is employed to predict the onset of buoyancy-driven motion. Under the principle of exchange of stabilities, the stability equations are derived in self-similar boundary-layer coordinate. The present predictions suggest the critical R D , and the onset time and corresponding wavenumber for a given R D . The onset time becomes smaller with increasing R D and follows the asymptotic relation derived in the infinite horizontal porous layer. Keywords Buoyancy-driven convection · Porous medium · Cylindrical geometry · Linear stability analysis

Nomenclature a∗ C c De g K P R (r, θ, Z ) (¯r , θ, z) RD

√ Modified wave number, αl,m τ Concentration (M) Dimensionless concentration Effective diffusivity (m2 /s) Gravitation acceleration (m/s2 ) Permeability of porous media (m2 ) Pressure (Pa) Radius of cylinder (m) Cylindrical coordinates Dimensionless cylindrical coordinates Darcy–Rayleigh number, R D = gβ K C R/(ε De ν)

M. C. Kim (B) Department of Chemical Engineering, Jeju National University, Jeju 690-756, Republic of Korea e-mail: [email protected]

123

396

M. C. Kim

t U W w

Time (s) Velocity vector (m/s) Vertical velocity (m/s) Dimensionless vertical velocity

Greek Letters αl,m β ε ζ μ ν τ ρ

Wave number Volumetric expansion coefficient (M−1 ) Porosity of porous media Dimensionless similarity variable, z/τ 1/2 Viscosity (Pa s) Kinematic viscosity (m2 /s) Dimensionless time Density (kg/m3 )

Subscripts 0 1 c

Basic quantities Perturbation quantities Critical conditions

Superscript *

Transformed quantities

1 Introduction The buoyancy-driven convection is well-known natural phenomena, and has attracted many researchers’ interests due to a wide variety of engineering applications, such as geothermal reservoirs, agricultural product storage system, packed-bed catalytic reactors, the pollutant transport in underground and the heat removal of nuclear power plants. Recently the phenomena of natural convection in porous media have been investigated in connection with the enhanced carbon dioxide dissolution into the saline water confined within the geologically stable formations (Ennis-King et al. 2005; Xu et al. 2006; Riaz et al. 2006; Hassanzadeh et al. 2006), and the novel-enhanced oil recovery (EOR) method where diffusion–convection may be employed as drive to force heavy, viscous oil to move (Rashidi et al. 2000). The dissolution process can be enhanced by the motion driven by the density gradient. Therefore, the onset condition of the buoyancy-driven motion is important to understand this process. The onset of buoyancy-driven motion in the horizontal porous layer begins with Horton– Rogers–Lapwood convection (Horton and Rogers 1945; Lapwood 1948). Later, Wooding (1962) suggested the critical condition for the long cylinder with unstable and linear density profile theoretically. However, most of engineering applications relating to rapid changing of temperature and concentration fields involve developing and nonlinear basic field. Therefore it is important to predict when the buoyancy-driven motion sets in. Recently, Riaz et al. (2006) analysed the onset of convection in porous media under the time-dependent concentration

123

Theoretical Analysis of Convective Instability

397

Fig. 1 Schematic diagram of system considered here

Cu g

Z=0 initial sharp concentration profile

homogeneous porous medium 2R Cl

field by using the self-similar coordinate. They showed that stability analysis in the selfsimilar, boundary layer coordinates yields the accurate stability criteria for the infinite depth system. In this study, the stability of an initially quiescent, liquids-saturated, cylindrical long porous layer experiencing mass diffusion is considered. The onset of convective motion by the buoyancy forces is analyzed under the linear stability theory. New stability limits, such as R D,c and tc will be suggested.

2 Theoretical Analysis 2.1 Governing Equations The system considered here is an initially quiescent, cylindrical porous column saturated with two miscible liquids, as shown in Fig. 1. The homogeneous and isotropic porous medium is confined within the long cylinder of radius R. At the initial time t = 0, it is assumed that the interface is sharp and that the liquids have uniform solute concentration Cl (lower fluid) and Cu (upper fluid). As discussed by Wooding (1962), this initially sharp interface model may not be realistic, but it describes the experimental situation approximately. Under the Boussinesq approximation, the governing equations can then be written as follows (Riaz et al. 2006): ∇ · U = 0, μ U = −∇ P − ρβgC, K

ε

∂C + U · ∇C = ε De ∇ 2 C, ∂t

(2.1) (2.2) (2.3)

where U (= er U + eθ V + ez W ) is the Darcy velocity vector in the cylindrical (r, θ, Z )coordinate system, μ the liquid viscosity, K the permeability, P the pressure, g the gravitational acceleration, ε the porosity, C the solute concentration, t time, β (= (1/ρ) (∂ρ/∂C)) the volumetric expansion coefficient and De the effective molecular diffusivity in the porous medium. The boundary conditions for the velocity and concentration fields are

123

398

M. C. Kim

U · n = 0 and (UC − De ∇C) · n = 0 at r = R, U→0 and C → Cl as Z → −∞, U → 0 and C → Cu as Z → ∞,

(2.4)

where n represents normal vector. The boundary conditions represent no radial mass flux through the cylinder wall and no through-flow and fixed concentration at the lower and upper boundaries. For the system of which density profile is linear and time-independent, i.e. (∂ρ/∂ Z ) = constant, Wooding (1959) suggested the following critical condition of onset motion:   gK R 2 ∂ρ = 3.390. (2.5) ε De μ ∂ Z For the time evolving system, however, the motion can be onset before the density profile becomes linear. Therefore, the stability problem becomes time dependent and very difficult, and the critical time tc and the minimum number to mark the onset of buoyancy-driven motion have practical importance. For this transient stability analysis we define a set of nondimensionalized variables τ, z, c0 by using the scale of time R 2 /De , length R and concentration C = (Cu − Cl ). Then the basic conduction state is represented in dimensionless form of ∂c0 ∂ 2 c0 , = ∂τ ∂z 2 with the following initial and boundary conditions,  1 z c0 (0, z) = 0 z  1 z c0 (τ, z) → 0 z

(2.6)

>0 , ≤0

(2.7a)

→∞ → −∞

(2.7b)

By using a similarity transformation, the above equations can be solved as follows:    ζ 1 c0 = 1 + erf , 2 2 √ where ζ = z/ τ .

(2.8)

2.2 Stability Equations Under the linear stability theory, disturbance equations can be formulated, in dimensionless form, in term of the concentration component c1 and the vertical velocity component w1 by decomposing Eqs. (2.1–2.3):

where ∇ 2 = by

∂2 ∂z 2

∇ 2 w1 = −∇12 c1 , (2.9) ∂c0 ∂c1 (2.10) + R D w1 = ∇ 2 c1 , ∂τ ∂z   ∂2 + ∇12 , ∇12 = r1¯ ∂∂r¯ r¯ ∂∂r¯ + r¯12 ∂θ 2 . The proper boundary conditions are given w1 → 0 and c1 → 0 as z → ±∞, ∂w1 ∂c1 = = 0 at r¯ = 1, ∂ r¯ ∂ r¯

123

(2.11a) (2.11b)

Theoretical Analysis of Convective Instability Table 1 Some lowest zeros of αl,m from Eq. (2.14c) corresponding to the least stable modes

(l, m) 1,1 αl,m

399

2,1

0,1

3,1

4,1

1,2

5,1

2,2

0,2

1.841 3.054 3.832 4.201 5.318 5.331 6.416 6.706 7.016

The velocity component has the scale of De /R and the concentration component has the scale of ε De ν/(gβ K R). Here, R D is the Darcy–Rayleigh number based on the radius of the cylinder defined by RD =

gβ K C R , ε De ν

(2.12)

where ν denotes the kinematic viscosity of liquids. Since the equations are linear with respect to r¯ and θ , the following normal mode analysis can be applied (Chandrasekhar 1961): [w1 , c1 ] = [w¯ 1 (τ, z) , c¯1 (τ, z)] X lm (¯r , θ ) .

(2.13)

The cylindrical harmonics X lm satisfies the following relation(Wooding 1959):

where

2 ∇12 X lm = −αl,m X lm ,

(2.14a)

  X lm (¯r , θ ) = Jm αl,m r¯ exp (imθ ) , m = 0, 1, 2, . . . ,

(2.14b)

  Jm αl,m = 0, l = 1, 2, . . . ,

(2.14c)

and

where Jm is the first kind Bessel function. Several smallest zeros of Eq. (2.14c) are summarized in Table 1. The non-axisymmetric case of m > 0 satisfies the continuity requirement through the factor exp (imθ ) and the axisymmetric one of m = 0 also does through the integral over the cylinder cross-section as 1 0

    J1 αl,0 r J0 αl,0 r dr = = 0, αl,0

by the Eq. (2.14c) (Wooding 1959). Therefore, the stability equations are summarized as:  2  ∂ 2 2 − α c¯1 , ¯ 1 = −αl,m l,m w ∂z 2  2  ∂ c¯1 ∂c0 ∂ 2 + R D w¯ 1 = − α l,m c¯1 . ∂τ ∂z ∂z 2

(2.15)

(2.16) (2.17)

under the following boundary conditions, w¯ 1 → 0 and c¯1 → 0 as z → ±∞,

(2.18)

Our goal is to find the critical time τc for a given R D and critical R D,c to mark the onset of convection for by using Eqs. (2.16–2.18). ∗ For the infinite-porous  √ layer, it is natural that c¯1 (τ, z) = c (τ, ζ ) and w¯ 1 (τ, z) = w ∗ (τ, ζ ), here ζ = z/ τ the similarity variable introduced already in the base state of

123

400

M. C. Kim

Eq. (2.8). Using the relations of ∂ c¯1 /∂τ = ∂c∗ /∂τ −ζ /(2τ )∂c∗ /∂ζ , the base state and the perturbation equations can be expressed as:  2  ∂ ∗2 −a (2.19) w ∗ = a ∗2 c∗ , ∂ζ 2   ∂2 dc0 ∂c∗ ζ ∂ τ − + − a ∗2 c∗ = −R ∗D w ∗ , (2.20) ∂τ ∂ζ 2 2 ∂ζ dζ   dc0 1 (2.21) = √ exp −ζ 2 /4 , dζ 2 π with the following boundary conditions: w ∗ → 0 and c∗ → 0 as ζ → ±∞. √ √ Here a ∗ = αl,m τ and R ∗D = R D τ .

(2.22)

3 Solution Procedure 3.1 Spectral Analysis Using the generalized Fourier series, c1 can be expanded as c∗ =

∞

Ai (τ ) κi φi (ζ ),

(3.1)

i=0

where κi is the normalization factor discussed later, and φi (ζ ) is the eigenfunction of the operator L and they should satisfy the following Sturm–Liouville equation: Lφi = −λi φi ,

(3.2)

under the following boundary conditions, φi (−∞) = φi (∞) = 0.

(3.3)

The solution of the above Sturm–Liouville boundary value problem, Eqs. (3.2) and (3.3), is  2   ζ ζ φi = (−1)i Hi exp − (3.4a) 2 4 and i +1 , i = 0, 1, 2, · · · , (3.4b) 2

  i    is the ith Hermite polynomial. The above where Hi (ξ ) = (−1)i exp ξ 2 ddξ i exp −ξ 2 eigenfunctions are orthogonal, since λi =

∞



 ζ2 κi κ j φi φ j exp dζ = δi j , 4

(3.5)

0

   √ −1 where κi = is the normalization factor and exp ζ 2 /4 is the  (i + 1) 2n+1 π weighting function of the Sturm–Liouville equation (3.2).

123

Theoretical Analysis of Convective Instability

401

From Eqs. (2.19) and (3.1), w1 is expressed as w1 =

∞

  Ai (τ )κi ψi a ∗ , ζ ,

(3.6)

i=0

where ψi can be obtained by solving   2 d ∗2 ψi (ζ ) = a ∗2 φi (ζ ) , −a dζ 2

(3.7)

under the following boundary conditions: ψi → 0 as ζ → ±∞

(3.8)

Using a method of the variation of parameters (Wessel-Berg 2009) and an inverse operator technique (Kim and Choi 2012), Eqs. (3.7) and (3.8) can be solved as: ⎧ ⎫ ⎡ ⎧ ζ ⎫⎤ ⎪ ⎪ ∞  ⎨ ⎬ ⎨ ⎬ ∗         a ⎢ ⎥ exp −a ∗ ξ φi dξ +exp −a ∗ ζ exp a ∗ ξ φi dξ ⎦ .(3.9) ψi = − ⎣exp a ∗ ζ ⎪ ⎪ ⎩ ⎭ 2 ⎩ ⎭ −∞

ζ

This solution can be expressed recursively as ψi = 4a ∗2 (ψi−2 + φi−2 ) ,

i = 2, 3, 4 . . . ,

(3.10a)

with

       ∗   ∗2   ∗  a∗ √ ζ ζ ∗ ∗ ψ0 = − +a +exp −a ζ 1 + erf −a π exp a exp a ζ erfc 2 2 2 (3.10b)

and

! "  2 # $   ∗2   ∗  ζ ζ ∗ ∗√ ∗ +a +a ψ1 = a exp a exp a ζ exp − − a πerfc 2 2 ! "    2 # $%  ∗  ζ ζ ∗ ∗√ ∗ + exp −a ζ exp − −a −a − a π 1 + erf 2 2 ∗

(3.10c) Substituting the above solutions (3.1) and (3.6) into Eq. (2.20), we obtain  2 ∞ ∞ ∞   1 d Ai ζ λi + a ∗2 Ai κi φi − R ∗D √ exp − τ Ai κi ψi , (3.11) κi φi = − dτ 4 2 π i=0 i=0 i=0   By using the orthonormal property of φi (ζ ), multiplying Eq. (3.10) by exp ζ 2 /4 κi φi (ζ ) and integrating over the range of ζ , the following system of homogeneous linear equation for the constants Ai is obtained: dAi τ dτ



∞ κi2 φi2 exp −∞

 2  ∞   ζ2 ζ ∗2 κi2 φi2 exp dζ = Ai λi − a dζ 4 4 −∞

∞

R∗ − √D κi Ak κk 2 π k=0

∞ φi ψk dζ,

(3.12)

−∞

123

402

M. C. Kim

or in matrix form τ

da = Ba, dτ

  where Bi j = − λi−1 + a ∗2 δi−1, j−1 −

(3.13)

R ∗D √ κ κ h , hi j 2 π i−1 j−1 i j

=

&∞ −∞

φi−1 (ζ ) ψ j−1 (ζ ) dζ ,

and a = [A0 , A1 , A2 , . . .]T . After an integration by parts, the following relation is obtained: ∞ hi j =

1 φi−1 (ζ ) ψ j−1 (ζ ) dζ = − ∗2 a

0

∞



 D 2 − a ∗2 ψi−1 (ζ ) ψ j−1 (ζ )dζ

0

1 = − ∗2 a

∞

  1 ψi−1 (ζ ) D 2 −a ∗2 ψ j−1 (ζ )dζ = − ∗2 a

0

∞

  ψi−1 (ζ ) D 2 − a ∗2 ψ j−1 (ζ ) dζ

0

∞ φ j−1 (ζ ) ψi−1 (ζ ) dζ = h ji .

=

(3.14)

0

This means the matrix B is symmetric and therefore it is normal. Usually, in the functional analysis, the norm of the disturbance c1 is written as ⎡∞ ⎤1/2  2  ' ζ c1 = ⎣ c12 exp dζ ⎦ = a T a, (3.15) 4 0

where a = [A0 , A1 , A2 , . . .]. Recently, Kim and Choi (2012) showed that the definition of the norm does little effect on the stability criterion. To trace the growth of disturbance, we defined the growth rate σ ∗ as   1 d c1 1 T B + BT ∗ σ = = T a a. (3.16) c1 dτ a a 2 Since the matrix B is symmetric, i.e. B = BT , the growth rate is reduced to σ ∗τ =

1 T a Ba. aT a

(3.17)

For the normal matrix B, the following Rayleigh quotient can be defined (Amunson 1966): R (B, v) =

v T Bv . vT v

(3.18)

And, the maximum value of R (B, v) reaches the maximum eigenvalue of B when v is the corresponding eigenvector, i.e. R (B, v) ≤ max {eig (B)} .

(3.19)

By using Eqs. (3.17) and (3.18), the growth rate of the most dangerous disturbances at a given time τ can be expressed as: σ ∗ τ = max {eig (B)} at τ = τ,

(3.20)

and the corresponding eigenvector is the most dangerous disturbances. In this study, the neutral stability condition is determined by setting σ ∗ = 0.

123

Theoretical Analysis of Convective Instability

403

3.2 Quasi-Steady-State Approximation (QSSA) Traditionally, time-dependent system has been analyzed by employing the quasi-steady-state approximation (QSSA) such as the frozen-time analysis. Under the QSSA, all the quantities are frozen at a certain τ , and the disturbance quantities c1 and w1 are assumed to have the following forms: )   ( (3.21) [c1 , w1 ] = c1∗ (ζ ) , w1∗ (ζ ) exp σ ∗ τ at a given τ. Then, the stability Eqs. (2.19–2.21) reduced as  2  D − a ∗2 w ∗ = a ∗2 c∗ ,    2 1 ζ ζ D 2 + D − a ∗2 − σ ∗ τ c∗ = −R ∗D w ∗ √ exp − , 2 4 2 π

(3.22) (3.23)

under the boundary conditions (2.22). Since w ∗ and c∗ should have the same symmetry relation in order to fit the governing equation and the boundary conditions, we transform the boundary conditions to Dwe∗ = Dce∗ = 0 at ζ = 0,

(3.24a)

wo∗ = co∗ = 0 at ζ = 0.

(3.24b)

or The first case means that w ∗ and c∗ are even functions and the second one does that w ∗ and c∗ are odd ones. Since the even mode is more unstable than the odd one, which will be shown later, we will find even mode solution which satisfies the boundary condition (3.24). The stability Eqs. (3.22–3.24) are solved by employing the outward shooting scheme (Kim 2012). In order to integrate these stability equations the proper values of we∗ and ce∗ at ζ = 0 are assumed for a given a ∗ and R ∗D . Since the stability equations and their boundary conditions are all homogeneous, the value of we∗ (0) can be assigned arbitrarily and the value of the growth rate σ is assumed. This procedure can be understood easily by taking into account of the characteristics of eigenvalue problems. After all the values at ζ = 0 are provided, this eigenvalue problem can be proceeded numerically. Integration is performed from ζ = 0 to an arbitrary upper with the fourth order Runge– Kutta–Gill method. If the guessed value of σ and ce∗ (0) are correct, we∗ and ce∗ will vanish at the upper boundary. To improve the initial guesses the Newton–Raphson iteration is used. For the present system, the position of the upper boundary is infinity. To consider the infinity boundary effect, the Shanks transformation is used. 3.3 Initial Value Problem Approach (IVPA) Through the above spectral analysis, the partial differential Eqs. (2.19–2.22) are reduced into the simultaneous ordinary differential equations (3.13), without spatial discretization. Here, we employ the initial value problem approach (IVPA) to obtain the solution of Eq. (3.13). To solve the initial value problem of Eq. (3.13), the proper initial condition should be considered. For the limiting case of τ → 0, Bi j = −λi δi j since w ∗ → 0 can be assumed ∗ from Eqs. (3.6) and (3.10),  since  a → 0. In this case, the most unstable disturbance is ∗ 2 c1 (τ, ζ ) = A0 (τ ) κ0 exp −ζ /4 and its growth rate is σ ∗ τ → −1/2 as τ → 0.

(3.25)

In the present IVPA, this most unstable initial disturbance as an initial value for Eq. (3.13).

123

404

M. C. Kim

growth rate, σ *

1.0

0.5

0.0

-0.5

Equation (4.2) 1st approximation 3rd approximation numerical (QSSA)

-1.0 -4 10

-3

10

-2

τ

-1

10

10

Fig. 2 Growth rates calculated from the various methods for R D = 50 and a = α0,1 2 α 1,1

RD =50

α 2,1

growth rate, σ

α 0,1 α 3,1

1

τc 0

-1 -4

10

-3

10

-2

τ

10

-1

10

Fig. 3 Growth rates obtained from the spectral analysis for the cases of R D = 50 and several smallest wavenumber

4 Results and Discussion For the specific case of R D = 50 and a = α0,1 , the growth rates based on the various approximations are compared in Fig. 2. The curves obtained by the first, and third approximations are compared with the numerical shooting solution. As shown in the figure, all curves for the limiting case of τ → 0 converge to the one by Eq. (3.32). With increasing τ , the third approximations start to deviate from the first approximation, and the numerical shooting solution supports the third approximation. Considering this figure, from now on, we adopt the third approximation to obtain the map of the temporal growth rate for the wide range of wavenumbers. The effect of the wavenumber on the growth rate is considered in Fig. 3, where the growth rates for the several smallest wavenumber are given. The neutral stability condition

123

Theoretical Analysis of Convective Instability Fig. 4 Neutral stability curves for the various approximations. (a) R ∗D vs. a ∗ . (b)R D τ 1/2 vs. a/R D from the third approximation

405

R*D

10 2

10 1

Equation (4.3)

10 0 0.0

1st approximation 3rd approximation QSSA, even mode QSSA, odd mode

(a) 0.5

1.0

1.5

2.2

a* 40

third approximation

RDτ1/2

30

20

0.1172 10

(b) 0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

a/RD Table 2 Comparison of critical conditions among the approximations

R ∗D,c ac∗

First approximation

Second approximation

Third approximation

Numerical shooting

5.65

5.65

5.58

5.58

0.42

0.42

0.43

0.43

determined by setting σ ∗ = 0 is compared in Fig. 4. Also, this figure shows that the even mode is more unstable that odd one. The critical conditions which correspond to the minimum point of each curve are given in Table 2.The region above the curve denotes the unstable state whereas that below the curve is stable one. For the limiting case of a ∗ → 0, from Eq. (3.10), it can be assumed that ψ0 >> ψ2n and ψ1 >> ψ2n+1 for n = 1, 2, 3, . . .. Therefore, the characteristic matrix B can be truncated as  B=

 B11 0 , 0 B22

(4.1)

123

406

M. C. Kim 1.0

growth rate, σ

*

0.5

0.0

-0.5

-1.0 -4 10

spectral method numerical shooting (QSSA) IVPA -3

10

-2

τ

10

-1

10

Fig. 5 Comparison of the grow rates calculated from the various methods for R D = 50 and a = α0,1 . The third approximations are used for the spectral method and the IVPA

    R∗ R∗ where B11 = − 21 + a ∗2 − 2√Dπ κ02 h 11 and B22 = − 1 + a ∗2 − 2√Dπ κ12 h 22 . Combining Eqs. (3.20) and (4.1), the growth rate can be approximated as   1 1 ∗ ∗2 σ τ →− (4.2) + R ∗D a ∗ as a ∗ → 0. +a 2 2 √ Since κ0−2 = 2 π and h 11 → −2a ∗ π as a ∗ → 0. Using this relation, the neutral R ∗D which can be obtained by setting σ ∗ = 0 is approximated as   1 + 2a ∗2 ∗ as a ∗ → 0. (4.3) RD → a∗ As shown in Fig. 4, for the limiting case of a ∗ → 0, Eq. (4.3) represents the neutral stability curves quite well. Now, consider the effect of the solution method on the stability characteristics. In Fig. 5, the IVPA solutions for the specific case of R D = 50 and a = α0,1 are compared with those from the spectral method and the numerical solution based on the QSSA formulation. It is interesting that the present growth rate is insensitive to the solution methods. This consistency is found in the similar system whose characteristic matrix B is normal (Kim and Choi 2012). For the limiting case of τ → 0, the present methods reproduce the initial growth rate given in Eq. (3.25). Unlike the horizontally unbounded system of infinite horizontal layer, the cylinder wall restricts the feasible wavenumber through Eq. (2.14c), that is the wavenumber has discrete values. For the several smallest wavenumbers, the neutral stability conditions for a given R D is summarized in Fig. 6. Based on this figure, the critical time and the critical wavenumber fora given  R D are represented in Fig. 7. As shown in this figure, non-axisymmetric mode of J1 α1,1 r¯ exp (iθ ) gives the most unstable mode and suggests the minimum R D to induce the convective motion. The present critical condition is R D,c =

123

0.1172 = 15.71, α1,1

(4.4)

Theoretical Analysis of Convective Instability

407

Fig. 6 Critical time from the third approximation for the smallest nine wavenumbers 10 0 asymtote α 1,1 10 -1

τc

15.71

α 2,1 α 0,1 α 3,1

10 -2

-1

ac ~0.078R D

10 -3 τc=(5.58/RD ) 10 -4 1 10

2

2

10

3

10

RD Fig. 7 Stability diagram showing the critical R D , τc and corresponding a

  2 of Wooding (1959), where which is quite different from gK R 2 /ε De μ (∂ρ/∂ Z ) = α1,1 (∂ρ/∂ Z ) is assumed to be constant for a long cylinder. The critical time decrease with increasing R D and if all the possible wavenumbers are superposed, it follows the asymptotic relation of τc = (5.58/R D )2 which is derived in the unbounded horizontal domain. In case of high R D the buoyancy-driven convection can be set in at small time and therefore the buoyancy-driven instability is confined the narrow region near the phase interface. Therefore, the radial boundary effects can be neglected, and the system can be considered the infinite-horizontal layer. For the unbounded infinite-horizontal layer, the minimum value of R ∗D and corresponding a ∗ in the Fig. 5 can be used to determine the critical condition marking the onset of buoyancy-driven convection. The stability criteria can be expressed for the infinite-horizontal system based on the third approximation as follows:

123

408

M. C. Kim

R ∗D,c = 5.58 and a ∗ = 0.43.

(4.5)

Based on this result, for the large R D the onset time tc and corresponding wavelength λc can be expressed based on the propagation theory as follows: " #   1/2 2 εν De εν De for R D > 100. λc = 81.54 (4.6) tc = 31.14 gβ K C gβ K C

5 Conclusions The onset of buoyancy-driven motion in a liquid-saturated, cylindrical porous layer experiencing mass diffusion has been analyzed theoretically by using the linear stability theory. The linear stability theory predicts that the onset time of buoyancy-driven motion is a function of the Darcy–Rayleigh number based on the radius of cylinder. The radial boundary restricts the feasible wavenumber and makes an important role in the stability characteristics especially for a small R D . Based on the presetnt analysis, the minimum Darcy–Rayleigh number to onset the buoyancy-driven motion is determined as R D,c = 15.71 and the critical time is τc = (5.58/R D )2 for R D >> R D,c . Acknowledgments This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012R1A1A2038983).

References Amunson, N.R., Mathematical Methods in Chemical Engineering, Prentice Hall, Englewood Cliffs (1966) Chandrasekhar, S.: Hydrodynamic and Hydromagnetic Stability. Oxford University Press, Oxford (1961) Ennis-King, J., Preston, I., Paterson, L.: Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids 17, 084107 (2005) Hassanzadeh, H., Pooladi-Darvish, N., Keith, D.W.: Stability of fluid in a horizontal saturated porous layer: effect of non-linear concentration profile, initial, and boundary conditions. Trans. Porous Med. 65, 193–211 (2006) Horton, C.W., Rogers, F.T.: Convection currents in a porous medium. J. Appl. Phys. 6, 367–370 (1945) Kim, M.C.: Onset of radial viscous fingering in a Hele–Shaw cell. Korean J. Chem. Eng. 29, 1688–1694 (2012) Kim, M.C., Choi, C.K.: 2012. Linear stability analysis on the onset of buoyancy-driven convection in liquidsaturated porous medium. Phys. Fluids 24, 044102 (2012) Lapwood, E.R.: Convection of a fluid in a porous medium. Proc. Cambridge Phil. Soc. 44, 508–521 (1948) Rashidi, F., Bahrami, A., Soroush, H.: 2000 Prediction of time required for onset of convection in a porous medium saturated with oil and a layer of gas underlying the oil. J. Pet. Sci. Eng. 26, 311–317 (2000) Riaz, A., Hesse, M., Tchelepi, H.A., Orr Jr, F.M.: Onset of convection in a gravitationally unstable, diffusive boundary layer in porous media. J. Fluid Mech. 548, 87–111 (2006) Wessel-Berg, D.: On a linear stability problem related to underground CO2 storage. SIAM J. Appl. Math. 70, 1219–1238 (2009) Wooding, R.A.: The stability of viscous liquid in a vertical tube containing porous material. Proc. Roy. Soc. Lond. A 252, 120–134 (1959) Wooding, R.A.: The stability of an interface between miscible fluids. ZAMP 13, 255–266 (1962) Xu, X., Chen, S., Zhang, Z.: Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 29, 397–497 (2006)

123