The onset of Darcy-Bénard instability in a horizontal porous channel with a free surface using a thermal nonequilibrium model

The intent beyond this present work is to investigate the onset of instability in homogeneous horizontal porous channel saturated with Newtonian fluid bounded by free surface from the top and impermeable wall from the bottom. The upper free surface is cooled by a third kind boundary condition, while the lower one is heated by two different ways. The wall is modelled with uniform temperature in the first case and with a perfect heat flux (Model A) in the second case. Relaxing the assumption of local thermal equilibrium model can provide in general two different Biot numbers Bs and B f as well as two temperature profiles at the basic state, one for each phase. The linear stability of the basic motionless state is studied by using the normal modes analysis. The neutral stability curves are obtained numerically by applying a sixth order Runge-Kutta method in combination with a shooting method. The critical values of R at which the instability initiates in the porous medium with isothermal wall or isoflux wall are calculated.


Introduction
The Darcy-Bénard problem is the equivalent of Rayleigh-Bénard problem with a porous material.The classical version of Darcy-Bénard configuration is defined as a plane porous layer saturated with Newtonian fluid and confined between two parallel isothermal walls heated from below [1].The existence of the porous material in thermoconvective region implies a new expression of the critical Rayleigh number, which is reformulated as a product between the traditional Rayleigh number for the pure fluid and Darcy number.It is well-known as Darcy-Rayleigh number, R. By presence of the local thermal equilibrium condition, one energy balance equation is considered.In this case, the threshold value of the Darcy-Bénard instability is R = 4π 2 [4].The thermal boundary conditions can change to various combinations between Dirichlet, Neumann and Robin boundary conditions.The same can also happen for the case of the hydrodynamic boundary conditions.Following Nield and Bejan [4], the critical value R = 0, can be observed while the boundary conditions are subjected to a uniform pressure and heat flux.Moreover, changing the boundary conditions into isoflux or free surface forces the basic flow to be more unstable.On the other hand, the local thermal equilibrium can be relaxed when the temperature of the solid and fluid phases are not the same.In this case the energy equation for the porous medium is replaced by two equations, one for the solid phase and one for the fluid phase.Under this model, the criterion for the onset of convection in porous system depends on the effect of two nondimensional parameters: the inter-phase heat transfer parameter H, the thermal conductivity ratio of the fluid and solid phase γ.Broadly speaking, the local thermal equilibrium model can appear inside a porous medium either by a e-mail: lagziri-hajar7@hotmail.fr [2].Thus, to study the effects of isoflux boundary instead of isothermal with local thermal non equilibrium condition on the onset of instability, Yang and Vafai [3] modelled a standard formula of the heat flux by two mathematical approaches, called model A and model B. The first one defines the total heat flux as sum between Fourier heat flux of the solid and fluid phase, weighted by the volumetric fraction of each phase, while the second one equals the Fourier heat flux considered in the solid phase and in the fluid phase.Usually, model A is used for the case of impermeable wall with higher thermal conductivity [3].
The aim of this current work is to study the comparative behaviour between the effect of the two Biot numbers B s and B f provided by Robin boundary condition, on the threshold values of the onset of instability in the porous channel bounded by isothermal wall from the bottom as first case, then by isoflux wall as second case when the condition of the local thermal equilibrium is abolished.

Mathematical model
The setup of the problem studied here is sketched in Fig 1, defined as an horizontal saturated porous layer confined between a free surface with uniform pressure at the upper boundary and rigid impermeable wall at the lower one.The free surface is subjected to Robin boundary condition, in meantime the rigid wall is modeled as isothermal T w in the first case then as isoflux q w (Model A) in the second case.With the existance of the local thermal non equilibrium model, the difference between the thermal conductivity coefficient of the solid and fluid phase, and the coefficient of the heat transfer between the solid phase h s or the fluid phase h f with an external environment yield two Biot numbers, B s and B f as well as two energy balance equations coupled by inter-phase heat transfer coefficient h.The saturated porous layer is gouverned by Darcy's law and Boussinesq approximation.The curl operator is used in Darcy's law equation in order to neglect the pressure field.On the other hand, the governing equations are invariant by rotation around the z axes.Thus, the velocity field can be defined in plan (x, z) with the streamfunctions w = −ψ ,x and u = ψ ,z .Hence, the system is reformulated as The subscript "s" and " f " denote a solid and fluid phases while ", z" and ", t" refer to first derivative with respect to z axes and time.The dimensionless boundary conditions are The Eqs. ( 1) and ( 2) are thus non-dimensionalized by the following transformations ( where the new dimensionless variables are Here, α m is the effective thermal diffusivity, k m is the effective conductivity in porous medium, K is the permeability, L is the thickness of the porous layer, φ is the porosity of the medium, α s and α f are the thermal diffusivities of the solid and fluid phases.

A basic non equilibrium state
Broadly speaking, when the local thermal non equilibrium is applied, the distribution of the temperature profile of the fluid phase and solid skeleton in both cases are different everywhere in the porous medium, see for instance [5].Thus, to prescribe the temperature and velocity fields at the basic state, we solve Eqs.(1) by using Eqs.
(2) then we report that The temperature field for the isothermal case The temperature fields for the isoflux case The subscript "b" denotes the basic flow and Ω = H(1 + γ).
The dimensionless parameters H and γ presente the interphase heat transfer coefficient, and thermal conductivity ratio of the fluid and solid phases, while R is the Darcy-Rayleigh number.

Stability analysis
The basic state is perturbed with an arbitrarily small disturbance reformulated in the following way The parameter is a small perturbation, while θ f , θs and ψ are the dimensionless fluid phase temperature disturbance, the dimensionless solid phase temperature disturbance and the dimensionless streamfuction disturbance, respectively.In order to obtain linear equations we introduce the decomposition of Eqs.(10) into Eqs.( 1) and (2).Then, by taking into account << 1, we remove all the additional terms of O( 2) Furthermore, according to the normal modes method, the solution of the linearized equations are expressed as superposition of the plane waves, namly ψ(x, z, t) = i ψ(z)e i(ax−ωt) , θs (x, z, t) = θ s (z)e i(ax−ωt) , Such that ω = ω r + iω i .The complex parameter of ω is the growth rate while the real one is the angular frequency.The parameter a is the dimensionless wave number and ψ, θ s , θ f are the dimensionless disturbance amplitudes.Thus, by substituting Eqs.(12) into Eqs.(11), the ordinary differential equations can be obtained as with the relative boundary conditions Analytically, the variable coefficients defined in the temperature field of the fluid phase for both cases, Eqs. ( 6) and ( 9), prevent us to emphasize the validity of the principle of exchange of stabilities.
In the general case, specifically when the two Biot numbers are different, the principle of exchange of stabilities is proved only numerically by the absence of the oscillating branch in the thermoconvective instability which means ω r = 0. On the other hand, we may assume ω i = 0 since the system involves modes that describe the threshold condition of neutral stability, see [5].Thus, the ordinary differential equations with ω = 0 become where the subscript ", zz" refers to second derivative with respect to z axes.

Numerical solutions
The eigenvalue problem, Eqs.(15), described for the full case can be solved by using two numerical procedures based on sixth order Runge Kutta solver coupled with shooting method.However, to apply the solver function employed in the setup of the Runge Kutta method, we must add three addtional conditions at boundary z = 0 in order to complete the initial boundary conditions defined in Eqs. ( 14) The homogeneity displayed in the eigenvalue problem yields the condition ψ ,z (0) = 1.The unknown parameters c 1 and c 2 are assessed in the setup used by a shooting method, with the boundary condition defined at free surface.The numerical procedures employed here are implemented in Mathematica software (© Wolfram Research).The three last conditions prescribed at the target z = 1 allowing us to compute the parameters c 1 , c 2 and the eigenvalue R, for any assigned input values of (a, γ, H, B f , B s ).Moreover, the curves R(a) drawn for each value given to the parameters (a, γ, H, B f , B s ) yield the neutral stability curves.

Discussion
The numerical results displayed in Figs 2-4 are reported for special general case based on some specific conditions suggested by the authors, in order to optimize the number of the governing parameters data (γ, H, B f , B s ) used in the study of the thermoconvective instability for both cases.
On the other hand, giving more accurate details about the behaviour of each parameter data inside a porous layer.Thus, the porosity of the porous layer in this work is assumed as 1/2 while both external heat transfer coefficients at the upper free surface are supposed to coincide.These assumptions provide a convenient relationship between the two different Biot numbers and γ in the form of B s = γB f .Figure 2 displays a neutral stability curves for fixed value of H = 10, with different values of γ: The graphic on the left hand side is for the case of isothermal wall and the graphic on the right hand side is for the case of isoflux wall.It is worth to mention that, in the case of isoflux wall the neutral stability curves are more sensitive to the effect of the parameter γ than the isothermal wall.The Figure 3 exhibits the trends of R c versus B f drawn for various values of γ with H = 0.1, while Figure 4 presents the plots of R c versus B s , drawn for different values of γ and with H = 0.1.The dashed horizontal lines describe the asymptotic values defined for B s → ∞ and B f → ∞.The values of R c in the case of isoflux wall increase consequently as one of the Biot number increases and γ decreases, on the other hand, the values of R c for isothermal wall display a non monotonic decrease with the Biot numbers and γ.In general, we can infer that, in the configurations of B s or B f → 0, the system becomes more stable with the case of isothermal wall and unstable with the case of isoflux wall.

Conclusion
The setup of Darcy-Rayleigh problem modified with free surface from the top, cooled by third boundary conditions and heated by uniform temperature from the bottom as first case, then by uniform heat flux (model A) as second case are numerically investigated.The effect of the local thermal non equilibrium model in the porous layer yields two different dimensionless Biot numbers, one for the fluid phase and the other for the solid skeleton.The analysis of the thermoconvective instability is carried out for special general case defined as B s = γB f in both cases.The numerical results obtained by employing a shooting method in conjunction with a sixth order Runge-Kutta method are characterized by the following behaviour: -The dimensionless parameters B f and B s affect signficantly the temperature gradient of the medium in the two cases.-The basic flow is more stable in the first case than in the second one.

1 Figure 2 .Figure 3 .
Figure 2. General case B s = γB f : The neutral stability curves for different values of γ, with fixed H = 10 for the cases of isothermal wall (left frame) and isoflux wall (right frame)

Figure 4 .
Figure 4. General case B s = γB f : Plots of R c versus B f for assigned values of γ, with fixed H = 0.1. )