Nonhomogeneous porosity and thermal diffusivity effects on stability and instability of double-diffusive convection in a porous medium layer: Brinkman Model

1 Introduction

Double diffusive convection in porous media has been a focus for researchers, since it has applications in a number of areas, including geophysics, the enhanced recovery of petroleum reservoirs, the underground diffusion of chemical wastes, seabed hydrodynamics, and crystal growth, for more details see [1], [2] and the references therein.

The critical Rayleigh number of the natural convection in anisotropic porous layers was analyzed by [3] and [4] where there was anisotropy in permeability and anisotropy both in permeability and thermal diffusivity, respectively. [5] proposed a general method of treating this kind of problem, and they reported results of numerical calculation of natural convection heat transfer near the critical Rayleigh number, using the Galerkin method. They also compared the calculated Nusselt number with the experimental values. All of these studies cited above assumed Darcy’s law. One of the recent and important study is introduced by [6]. In this study, the critical Rayleigh number at the onset of natural convection was studied by linear stability analysis for high porosity anisotropic horizontal porous layers of uniformly arraying vertical thin circular wires stretched across a hot and a cold surface. Navier-Stokes equations including flow resistance by the wires were solved since Darcy’s law cannot be applied due to high porosity.

Capone et al. [7] analyzed the onset of convection in a horizontal fluid-saturated porous layer uniformly heated from below in which the permeability and thermal diffusivity vary linearly or exponentially in the vertical direction. They studied this problem with constant gravity. However, Harfash [8] studied this problem with variable gravity with respect to the vertical direction and tested the validity of both the linear instability and global nonlinear energy stability thresholds for this problem using a three-dimensional simulation. Moreover, a model for double-diffusive convection in anisotropic and inhomogeneous porous media using Darcy’s law has been analysed in [9]. The effects of variable permeability, thermal diffusivity and variable gravity with respect to the vertical direction, have been studied in [9].

There is increasing interest in micro-electro-mechanical-systems (MEMS), and flow in microfluidic channels due to their applications in the electronics and related industries. In particular, at nanoscales there is increasing evidence that boundary conditions of slip type are needed rather than those of no-slip, cf. Cercignani [10], Duan [11], Duan & Muzychka [12], Lauga et al. [13], Morini et al. [14], Priezjev [15], Zhang et al. [16], Zhang et al. [17]. An especially important application of microscale flow involving slip boundary conditions is to flow in porous metallic foams. Lefebvre et al. [18] give many industrial examples of this and provide a thorough review of the state of the art. The object of this paper is to investigate a double diffusive convection problem in an inhomogeneous Brinkman porous media. In particular, the effect of variable permeability, thermal diffusivity and slip boundary conditions, has been studied.

Microchannels heat sinks with pin fins or ribs are considered as porous media. If the distance between the solid materials in micro dimensions, the slip is obtained in the boundaries (Knudsen number is appareled) for the flow of gases in the pores. The Knudsen number (Kn) is a dimensionless number defined as the ratio of the molecular mean free path length to a representative physical length scale. This length scale could be, for example, the radius of a body in a fluid. Knudsen number is expressed by the following equation.

where Kn Knudsen number, λ mean free path (m), D_h hydraulic diameter (m).

The mean free path of a particle, such as a molecule, is the average distance the particle travels between collisions with other particles. The mean free path is a function of temperature and pressure and is defined as follows

where K_b, T, P, and d are the Boltzmann constant 1.3806503E – 23 (m²kgs^–2K^–1), temperature (K), pressure (N/m²), and molecular diameter (m).

As Knudsen number increases, rarefaction effects become more important and continuum approach breaks down. However, it is convenient to differentiate the flow regimes in function of Kn, and the following classification:

1. For Kn < 10^–3, the flow is a continuum flow and it is accurately modeled by the compressible Navier-Stokes equations with classical no-slip boundary conditions.

2. For 10^–3 < Kn < 10^–1, the flow is a slip flow and the Navier-Stokes equations remain applicable, provided a velocity slip and a temperature jump are taken into account at the walls. These new boundary conditions point out that rarefaction effects become sensitive at the wall first.

3. For 10^–1 < Kn < 10, the flow is a transition flow and the continuum approach of the Navier-Stokes equations is no longer valid. However, the intermolecular collisions are not yet negligible.

4. For Kn > 10, the flow is a free molecular flow and the occurrence of intermolecular collisions is negligible compared with the one of collisions between the gas molecules and the walls.

In the gas reservoir (petroleum reservoirs), there are a difference between the pressures (temperature) of reservoir layers which makes the diffusion is transfer from high to low through the pores, where the most of pores are in the micro dimension or less which produces a slip in the boundaries.

The really interesting situation from both a geophysical and a mathematical viewpoint arises when the layer is simultaneously heated from below and salted from below. In this situation heating expands the fluid at the bottom of the layer and this in turn wants to rise thereby encouraging motion due to thermal convection. On the other hand, the heavier salt at the lower part of the layer has exactly the opposite effect and this acts to prevent motion through convective overturning. Thus, these two physical effects are competing against each other. Due to this competition, it means that the linear theory of instability does not always capture the physics of instability completely and (sub-critical) instabilities may arise before the linear threshold is reached. Due to the possibility of sub-critical instabilities occurring, it is very important to test the accuracy of the linear instability threshold in capturing the onset of the instability of the double diffusive problem when the layer is heated and salted from below. Recent contributions on the study of convective instabilities in fluid and porous media include [19,20,21,22,23,24,25,26,27,28]. In the last few years and via a new methodology which is called the auxiliary system method, Rionero proved the coincidence between linear and nonlinear stability thresholds for the convection problem in many real situations (for more details, see [29,30,31,32,34]).

The plan of the paper is as follows. In the next section, we present the governing equations of motion and derive the associated perturbation equations. Then, we perform linear instability analysis (Section 3) and global nonlinear stability analysis (Section 4) and derive numerically the instability/stability thresholds. Since the stability analyses involve eigenvalue problems with non-constant coefficients these problems must be solved numerically and a suitable numerical method is described in section 5. In Section 6, the numerical results for the linear theory and a direct comparison with those of the global nonlinear theories are presented.

3 Linear instability

In order to study linear instability, we discard the nonlinear terms in (2.11). A time dependence such as u_i = e^σtu_i(x), π = e^σtπ(x), θ = e^σtθ(x), ϕ = e^σtϕ(x) is now assumed and then, after removing, the pressure perturbation the linearized instability equations that arise from (2.11) are found to be

where Δ^* is the horizontal Laplacian Δ^* = ∂²/∂ x² + ∂²/∂ y², D = d/dz, and (3.1) hold on ℝ² × (0, 1). To proceed from equations (3.1), a plane tiling form h(x,y) is introduced, see e.g. [42, 43]), and then we put w = W(z) h(x,y), θ=Θ(z) h(x,y) and ϕ = Φ(z) h(x,y) and introduce the wavenumber a by Δ^*h = –a²h. Equations (3.1) then yield the eigenvalue problem

In (3.4), z ∈ (0, 1), and the boundary conditions are

Detailed numerical results are presented in Section 7.

4 Nonlinear energy stability theory

When adopting a linear analysis approach, the perturbation to the steady state is assumed to be small, and so nonlinear terms in the governing set of partial differential equations are discarded. It has been proved that linear analysis often provides little information on the behavior of the nonlinear system [42], so in such cases only instability can be deduced from the linear thresholds, as any potential growth in the nonlinear terms is not considered.

Let V be a period cell for a disturbance to (2.11), and let ‖⋅‖ and 〈⋅, ⋅ 〉 be the norm and inner product on L²(V). Next, multiply (2.11)₁, (2.11)₃ and (2.11)₄ by u_i, θ and ϕ and integrating over V. By introducing a coupling parameters λ₁, λ₂ > 0 we may then derive

where the functions E and I are given by

and

with the dissipation D being defined by

where dS is a surface element, Ψ₁(z)= 1+λ1g(z),Ψ₂ = 1 – λ₂ and u is explicitly written as u = (u, v, w). Define R_E by

where H is the space of admissible functions, i.e. u_i, θ, ϕ ∈ H¹(V) with u_i solenoidal and u_i, θ, ϕ satisfying the boundary conditions. Then from (4.1) we derive

Poincare’s inequality ensures that there is a constant c > 0 such that 𝓓 > cE and then if R_E > 1, one may show from inequality (4.6) that E decays exponentially and we have global nonlinear stability, i.e. for all initial data. The nonlinear stability threshold then requires the solution of (4.5). The Euler-Lagrange equations which arise from (4.5) are

where ζ is a Lagrange multiplier. To remove the Lagrange multiplier we take the third component of the double curl of (4.5)₁, and introducing the normal mode representation and notation as presented in section 3, thus (4.5) then becomes

together with boundary conditions (3.5). Numerical results are presented in Section 7.

5 Numerical technique

In this section, we use the Chebyshev collocation method to solve the eigenvalue system (3.4) and (4.8). In the Chebyshev collocation method, system (3.4) is rewritten in terms of second and third order derivatives only. Letting Π = DW, (3.4) can be expressed as the four 2nd order equations. The system is then transformed onto the Chebyshev domain (–1, 1) and the solutions W and θ treated as independent variables and expanded in a series of Chebyshev polynomials

then, we insert (5.1) into the equations (3.4), and then substitute the Gauss-Labatto points which are defined by

Thus, we obtain 4N – 4 algebraic equations for 4N + 4 unknowns W₀,..., W_N, Π₀,..., Π_N, Θ₀,..., Θ_N, Φ₀,..., Φ_N. Now, we can add six rows using the boundary conditions (3.5) as follows

The inner product of each equation is taken with some T_k and the orthogonality of the Chebyshev polynomials exploited to obtain the following generalised eigenvalue problem:

where X = (W₀, …, W_N, Θ₀, …, Θ_N, Φ₀, …, Φ_N), O is the zeros matrix, I(n₁, n₂) = T_n2(z_n1), D(n₁, n₂) = Tn2′ (z_n1), D²(n₁, n₂) = Tn2″ (z_n1), Y(n₁, n₂) = (1/g(z_n1))I(n₁, n₂), Σ(n₁, n₂) = f²(z_n1)I(n₁, n₂), Ω₀ = χa⁴f²(z_n1)I(n₁, n₂) + a²f(z_n1)I(n₁, n₂), Ω₁ = χf²(z_n1)(8D³(n₁, n₂) − 4a²D(n₁, n₂)) − 2f(z_n1)D(n₁, n₂) + f′(z_n1)I(n₁, n₂), Ω₂ = g(z_n1)(4D²(n₁, n₂) − a²I(n₁, n₂)) + 2 g′(z_n1)D(n₁, n₂), Ω₃ = 4D²(n₁, n₂) − a²I(n₁, n₂) n₁ = 0, …, N − 2 and n₂ = 0, …, N. We computed the differentiation matrices, which are corresponded to the trail functions (5.1) analytically using Matlab routines. For more details about the Chebyshev collocation method, see [44].

We have found the solution to (5.3) is achieved very accurately by employing the QZ algorithm of Moler and Stewart [45]. This algorithm relies on the fact that there are unitary matrices Q and Z such that QAZ and QBZ are both upper triangular. The algorithm then yields sets of values α_i, β_i which are the diagonal elements of QAZ and QBZ. The eigenvalues α_i of (5.3) are then obtained from the relation α_i = α_i/β_i, provided β_i ≠ 0. This is very important, since the way we have constructed B means it contains a singular band and we find one third of the β_i = 0; these β_i must be filtered out. Indeed, with the technique advocated here one ought always to consider the α_i and β_i, since as [45].

We have solved system (5.3) for eigenvalues σ_j by using the QZ algorithm from Matlab routines. Once the eigenvalues σ_j are found we use the secant method to locate where σjR,σj=σjR+σjI being the real and imaginary parts of eigenvalue σ_j. The value of R which makes σ1R=0,σ1R being the largest eigenvalue, is the critical value of R for a² fixed. We then use golden section search to minimize over a² and find the critical value of R² for linear instability. Numerical results are reported in the next section. In our use of the Chebyshev collocation method, we used polynomial of degree between 20 and 30. Usually 25 was found to be sufficient but convergence was checked by varying the degree by examining the convergence of the associated eigenvector (which yields the approximate associated eigenfunction).

Returning to the nonlinear eigenvalue system (4.8), the the Chebyshev collocation method yields:

where

Then, we can determine the critical Rayleigh Ra_E for fixed a², λ₁ and λ₂ . Next, we employ golden section search to minimize in a² and then maximize in λ₁ and λ₂ to determine Ra_E for nonlinear energy stability,

where for all R² < Ra_E we have stability. The derivatives in the optimization problem are not known and so a quasi-Newton (or similar) technique that does not require derivatives would have to be used. We have found the golden section search method works well on both the maximum and minimum problems. It may be a little more expensive in computer time than other techniques, but for problems like (5.5) it has been very reliable.

6 Stability analysis results

The numerical results are presented for g(z) = e^Az, g(z) = e^Bz, and N_L = N_U = N₀. The results in this paper are given for N_L = N_U = N₀, ε̃ =5, and ϑ = Le = 1. This selection of parameter values was corresponds to the values in previous studies which have been deal with the problem of nonhomogeneous porosity and thermal diffusivity effects on stability and instability of double-diffusive convection [7, 8].

The thresholds of both the numerical linear instability and nonlinear stability results are presented in Figure 1. More comparable linear and nonlinear thresholds are apparent as the onset of convection predicted by the linear theory becomes fully stationary. However, their agreement does deteriorate as the solute Rayleigh number becomes large, indicating that the linear theory may fail to suitably emulate the physics of the onset of convection. The behaviour of the linear instability curves is in good agreement with that seen in [46]. The kink in the curves represents the point at which convection switches from steady convection (σ = 0) to oscillatory (σ_r = 0, σ_i ≠ 0). Note that as R_c increases the onset of convection is more likely to be via oscillatory convection as opposed to steady convection. We observe from Figure 1 that for R_c = 0 there is very good agreement between the nonlinear stability and linear instability bounds. We also note that as R_c is small the linear and nonlinear results become closer and remark that since the critical Rayleigh numbers are so close (for R_c sufficiently small) we expect the linear analysis to have captured the essential physics of the onset of convection. However as R_c increases the agreement between the two thresholds is not so good and we highlight this area as a region of possible subcritical instabilities.

Fig. 1

Visual representation of stationary linear instability (solid line) and nonlinear stability (dashed line) thresholds, with critical Rayleigh number plotted against R_c, where N₀ = 1, χ = 0.1 and A = B = 1.

Figure 2 explores this concept in more detail by giving a visual representation of the linear instability thresholds for a variety of χ values for fixed A, B,ε, N₀ and R_c. Interesting, it is clear from figure 2 an increase in χ values causes the system to become slightly more stable. As the linear instability and nonlinear stability results clearly show slightly good agreement, we can conclude that (for the parameter ranges explored) the linear theory not accurately encapsulates the physics of the onset of convection. Thus, the novel instability results are supported by the considerable region of subcritical instabilities. However, for a fixed value of A = B = 1, N₀ = 1 and R_c = 15, we see that oscillatory convection is occurring for χ ≤ 0.575 and when χ > 0.575 we shall witness stationary convection. In Figure 3, we display the critical Rayleigh numbers at which stability and instability begin as a function of the slip coefficient, N₀. This Figure is given for A = B = 1, ε = 0.5, χ = 0.1 and R_c = 15. It is clear from Figure 3 that the region of potential subcritical instabilities between the linear instability and nonlinear stability thresholds is considerable. From Figure 3, we see that the oscillatory convection is always dominant in that the instability curve always lies below the stationary convection one for N₀ ≥ 0.05. Figure 4 shows how increasing A, corresponds, to destabilization. It is clear from Figure 4 that an increase in A causes the system the become more unstable, which we would physically expect. Again, it is very noticeable that the nonlinear energy stability curves are close to those of linear instability. This is reinforcing the fact that the linear curves are true representation that the physics of the onset of convection is being correctly reflected. The gap between the curves represents the small band where sub-critical bifurcation may possibly occur. For A < –1.4, the stationary convection become dominator in the linear instability thresholds. However, for A ≥ –1.4, the oscillatory convection appears in the linear instability.

Fig. 2

Visual representation of stationary linear instability (solid line) and nonlinear stability (dashed line) thresholds, with critical Rayleigh number plotted against χ, where A = B = 1, N₀ = 1 and R_c = 15.

Fig. 3

Visual representation of stationary linear instability (solid line) and nonlinear stability (dashed line) thresholds, with critical Rayleigh number plotted against N₀, where A = B = 1, χ = 0.1 and R_c = 15.

Fig. 4

Visual representation of stationary linear instability (solid line), oscillatory linear instability (dashed line) and nonlinear stability (dotted line) thresholds, with critical Rayleigh number plotted against A, where B = 1, ε = 0.5, N₀ = 1, χ = 0.1 and R_c = 14.

Figures 5 gives a visual representation of the linear instability and nonlinear stability thresholds, with critical thermal Rayleigh number Ra plotted against B. The remaining parameters are held fixed at A = 1, N₀ = 1, χ = 0.1 and R_c = 14. This figure shows the effect of increasing B on the critical Rayleigh number. It is clear from this Figure that an increase in B causes the system the become more stable, which we would physically expect. For B < –0.1, the stationary convection become dominator in the linear instability thresholds. When B ≥ –0.1, the oscillatory modes become present in the linear instability thresholds. Figure 5 demonstrates that Ra increases with increasing B which shows the stabilizing effect of B. It is very noteworthy that the nonlinear stability curves are close to those of linear theory. This shows that possible sub-critical instabilities may arise in a small range of Rayleigh numbers, and it also demonstrates that linear instability theory does correctly capturing the physics of the onset of convection.

Fig. 5

Visual representation of stationary linear instability (solid line), oscillatory linear instability (dashed line) and nonlinear stability (dotted line) thresholds, with critical Rayleigh number plotted against B, where A = 1, N₀ = 1, χ = 0.1 and R_c = 14.

7 Conclusions

In this paper we have explored the effect of slip boundary conditions on the double-diffusive convection in anisotropic and inhomogeneous porous media using the Brinkman model, utilising linear stability analysis and non-linear stability by means of energy functional. A comparison between the linear stability thresholds and energy stability thresholds is made. In both cases the thermal Rayleigh number is evaluated for different combinations of the flow governing parameters. The results indicate that the increasing in the slip coefficient N₀ has a destabilizing effect in the linear and nonlinear cases. Moreover, the results show that an increase in the permeability coefficient A causes a strong destabilization effect on the results. However, with increasing the thermal diffusivity coefficient B, the results show the stabilizing effect on the horizontal Rayleigh numbers in the linear and nonlinear cases.