Interaction of the Longwave and Finite-Wavelength Instability Modes of Convection in a Horizontal Fluid Layer Confined between Two Fluid-Saturated Porous Layers

The onset of convection in a three-layer system consisting of two fluid-saturated porous layers separated by a homogeneous fluid layer is studied. It is shown that both a longwave convective regime developing in the whole system and a finite-wavelength regime of convection concentrated in the homogeneous fluid layer are possible. Due to the hydraulic resistance of the porous matrix, the flow intensity in the longwave convective regime is much lower than that in the finite-wavelength regime. Moreover, it grows at a much slower pace with the increase of the Grashof number. Because of that, the long-wave convective regime becomes unstable at small supercriticalities and is replaced by a finite-wavelength regime.


Introduction
At the onset of convection in layers, the scale of convective cells is usually close to the transversal size of the channel.The original situation is discussed in Reference [1], where the onset of convection in a three-layer system consisting of two fluid-saturated porous layers separated by a pure fluid layer is studied.It was found that, there is a range of parameters where the neutral curves are bimodal, i.e., both, a longwave convection, taking place mainly in porous medium, and a finite-wavelength convection, concentrated mainly in a fluid with a characteristic horizontal scale close to the thickness of a layer, may coexist.Later on, the same phenomena were discovered and analyzed for two-layer systems made up of a superposed pure fluid layer and a fluid-saturated porous layer [2][3][4][5][6].In Reference [7], it is shown that high frequency transversal vibrations cause a stabilizing effect on the disturbances with any wave numbers, but affects finite-wavelength disturbances much more strongly than longwave ones.The present work is devoted to the investigation of finite-wavelength and longwave instability modes interaction in the supercritical parameter range for the three-layer configuration considered in Reference [1].The studies of the onset of convection in multi-layer systems consisting of one fluid layer and two fluid-saturated porous layers are important, due to their application to the directional solidification of binary alloys where a two-phase porous zone called the mushy zone is formed between the melt and the crystal.As shown in Reference [8], there are two types of instability in a two-layer system consisting of the melt zone and the mushy zone: short-wave instability and long-wave instability.

Governing Equations and Boundary Conditions
The problem of convection in a system which consists of an infinite horizontal layer of an incompressible viscous fluid with the thickness 2d embedded between two plane layers of porous medium saturated with the same fluid, each with a thickness of h m , is considered (Figure 1).The temperature of the external boundaries is fixed, and the temperature of the bottom boundary is higher than that of the top (heating from below).The origin of the coordinate system is set in the middle of the fluid layer.

Governing Equations and Boundary Conditions
The problem of convection in a system which consists of an infinite horizontal layer of an incompressible viscous fluid with the thickness 2d embedded between two plane layers of porous medium saturated with the same fluid, each with a thickness of m h , is considered (Figure 1).The temperature of the external boundaries is fixed, and the temperature of the bottom boundary is higher than that of the top (heating from below).The origin of the coordinate system is set in the middle of the fluid layer.Equations for thermal buoyancy convection in fluid in the Boussinesq approximation look like [9]: For the description of convective filtration in a porous medium, we use the Darcy-Boussinesq model [10]: Here, v  and u  are the convective flow velocity in fluid and the convective filtration velocity in a porous medium, respectively; T and ϑ are the temperature deviations in fluid and in the porous medium from a constant average temperature; m is the porosity coefficient; and K is the permeability coefficient of porous medium.Other designations are usual.The quantities marked by the index f concern the fluid, and those marked by the index m concern the porous medium.
Let us discuss the boundary conditions.The external boundaries are considered to be impenetrable, i.e., the vertical component of convective filtration velocity vanishes at these boundaries.At the same time the horizontal components of the filtration velocity on the external boundaries are generally non-zero, such that the filtration of fluid in the porous medium along the impenetrable boundary is possible.The temperature of the external boundaries is supposed to be fixed.Thus, the conditions at the external boundaries are as follows: Equations for thermal buoyancy convection in fluid in the Boussinesq approximation look like [9]: For the description of convective filtration in a porous medium, we use the Darcy-Boussinesq model [10]: Here, → v and → u are the convective flow velocity in fluid and the convective filtration velocity in a porous medium, respectively; T and ϑ are the temperature deviations in fluid and in the porous medium from a constant average temperature; m is the porosity coefficient; and K is the permeability coefficient of porous medium.Other designations are usual.The quantities marked by the index f concern the fluid, and those marked by the index m concern the porous medium.
Let us discuss the boundary conditions.The external boundaries are considered to be impenetrable, i.e., the vertical component of convective filtration velocity vanishes at these boundaries.At the same time the horizontal components of the filtration velocity on the external boundaries are generally non-zero, such that the filtration of fluid in the porous medium along the impenetrable boundary is possible.The temperature of the external boundaries is supposed to be fixed.Thus, the conditions at the external boundaries are as follows: Different types of boundary conditions at the interface of pure fluid and fluid-saturated porous medium were suggested (see Reference [10]).The conditions proposed by Beavers and Joseph [11] determine the jump of the tangential components of the velocity at the interface while the normal Fluids 2017, 2, 39 3 of 8 component of the velocity and pressure are assumed to be continuous.The conditions obtained by Ochoa-Tapia and Whitaker [12,13] by the direct averaging of the Navier-Stokes equations at the pore-scale level impose the continuity of flow velocities and normal stresses, while the tangential viscous stresses have a jump across the interface, describing the resistance of the porous matrix in the boundary layer, which has a thickness to of the order of K 1/2 .These conditions, as well as the Beavers-Joseph conditions, include an empirical parameter -the stress jump coefficient.The influence of the stress jump coefficient on the onset of thermal buoyancy convection in horizontal stratified fluid/porous layers was studied in Reference [14].It was shown that the increase of the stress jump coefficient strongly influences the fluid mode, inducing a more unstable situation at large wave numbers, whereas the porous mode remains unchanged.In Reference [15], the onset of thermal buoyancy convection in a system consisting of a fluid layer overlying a homogeneous porous medium was studied in the framework of the model including the Brinkman term, in a two-domain approach.A comparison of neutral curves with those obtained in the one-domain approach and in the framework of the Darcy model with the two-domain approach shows that the inclusion of the Brinkman term plays a secondary role in the stability results.
At zero value of this parameter, the tangential stress jump condition is reduced to the tangential stress continuity condition proposed in Reference [1].In References [16,17], the improvement of the approach by Ochoa-Tapia and Whitaker was suggested.
In our work, we use the conditions at the interface of pure fluid and fluid-saturated porous medium suggested by Lyubimov and Muratov in Reference [1], which include the continuity of temperature, heat flux, normal component of velocity and pressure, and the condition of vanishing of the tangential components of fluid velocity: The application of the boundary condition v x = v y = 0 is justified by the fact that the velocities of the convective filtration in a porous medium are generally small because of the small typical values of the porous medium permeability.In this case, the condition of the normal stress balance is reduced to pressure continuity.The advantages of the boundary conditions (Equation ( 4)) are their simplicity and the fact that they do not include any empirical parameters.
It is possible to show that, just as in the case of a homogeneous fluid, the condition for a conductive state of fluid in the described system is the linear dependence of temperature on vertical coordinate: Then, from the condition of heat flux continuity follows: and the condition: Further consideration is limited to the case when the thermal properties of the fluid and the "skeleton" are identical, that is: It is convenient to present the condition of pressure continuity on the interface of the porous medium and fluid in the terms of velocities.Using horizontal projections of momentum and continuity equations, we obtain this condition in the form: The analysis was restricted to the case of two-dimensional flows.In this case, it is convenient to introduce the stream function and vorticity for the description of the flow in the fluid layer and the stream function for the flow in the fluid-saturated porous layers.Equations rewritten in terms of stream function and vorticity in the dimensionless form are: where ∂x , and the last condition at z = ±1 is obtained from Equation ( 8), taking into account the smallness of the Darcy number.
The problem is characterized by four dimensionless parameters: the ratio of the porous layer thickness to the pure fluid layer thickness, the Darcy number, the Prandtl number, and the Grashof number:

Numerical Results
The problem was solved numerically by the finite difference method.The explicit finite difference scheme of the second order of approximation for spatial variables was used.The grid with the spatial step h = 1/20 was chosen for the main calculations, after the test calculations confirmed the convergence of the numerical results with the decrease of h.
As shown in Reference [1], the coexistence of finite-wavelength and longwave instability takes place at sufficiently large values of h in some range of Da values, which are extended while h increases.
In Figure 2, the neutral curves C(k) where C = √ Gr are plotted for h = 10 and different values of Da.It can be seen that at Da = 10 −3 and Da = 2 • 10 −3 the neutral curves are bimodal; moreover, at Da = 10 −3 , the finite-wavelength and longwave minima of the neutral curve correspond to nearly the same values of the Grashof number (the longwave minimum is just slightly lower).
On the basis of these results, the parameter values h = 10, Da = 10 −3 , Pr = 1 were chosen for the nonlinear calculations.The length of the computational domain was found to be equal to 10π, that is, approximately equal to the wavelength of the most dangerous longwave disturbances.
On the vertical boundaries of the area, the conditions of periodicity were imposed.Such conditions allow not only longwave solutions but also finite-wavelength ones, with the wavelength equal to 10π/n, where n is an integer value.
On the vertical boundaries of the area, the conditions of periodicity were imposed.Such conditions allow not only longwave solutions but also finite-wavelength ones, with the wavelength equal to 10 / n π , where n is an integer value.Naturally, it is easier to obtain finite-wavelength solutions in calculations with short computational domain.However, in this case we lose an opportunity to look after their interaction with longwave disturbances.As the symmetry of longwave and finite-wavelength solutions is the same, such an interaction may be sufficiently strong.
The numerical calculations show that, as expected, finite-wavelength stationary solutions describe the convective flow concentrated mainly in the fluid (Figure 3), and longwave solutions describe the flow spreading both through the layer of homogeneous fluid and the porous layers (Figure 4).As follows from the results of the calculations obtained in the framework of linear theory, critical values of the Grashof number for disturbances with n = 4, 5, 6, 7 are close to those for longwave disturbances with n = 1, although they are a little bit higher.
Naturally, it is easier to obtain finite-wavelength solutions in calculations with short computational domain.However, in this case we lose an opportunity to look after their interaction with longwave disturbances.As the symmetry of longwave and finite-wavelength solutions is the same, such an interaction may be sufficiently strong.
The numerical calculations show that, as expected, finite-wavelength stationary solutions describe the convective flow concentrated mainly in the fluid (Figure 3), and longwave solutions describe the flow spreading both through the layer of homogeneous fluid and the porous layers (Figure 4).
On the vertical boundaries of the area, the conditions of periodicity were imposed.Such conditions allow not only longwave solutions but also finite-wavelength ones, with the wavelength equal to 10 / n π , where n is an integer value.Naturally, it is easier to obtain finite-wavelength solutions in calculations with short computational domain.However, in this case we lose an opportunity to look after their interaction with longwave disturbances.As the symmetry of longwave and finite-wavelength solutions is the same, such an interaction may be sufficiently strong.
The numerical calculations show that, as expected, finite-wavelength stationary solutions describe the convective flow concentrated mainly in the fluid (Figure 3), and longwave solutions describe the flow spreading both through the layer of homogeneous fluid and the porous layers (Figure 4).The intensity of finite-wavelength solutions near the threshold grows according to the square root law, i.e., the excitation of short-wave convection occurs through supercritical bifurcation (Figure 5).
The calculations show that in the longwave mode, where the flow covers not only the fluid layer but also the porous medium, the flow intensity is lower than in the finite-wavelength mode, where the flow is concentrated in the fluid.This is caused by the hydraulic resistance of the porous medium skeleton.Moreover, the intensity of longwave solutions grows very slowly with the increase of the Grashof number.This is why already at small supercriticalities with respect to finite-wavelength disturbances, the intensity of the latter surpasses the intensity of the longwave solution (Figure 5).Apparently, the fast loss of stability of the longwave solution is related to this circumstance.Let us discuss the interaction of longwave and finite-wavelength disturbances at small supercriticalities.It turns out that the critical value of the Grashof number for the longwave instability mode is smaller that for the finite-wavelength mode.Because of this, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size: at the Grashof number values larger than 80 after the transient stage, depending on the initial conditions, a 4, 5, or 6 wave regime is established.
The process of the structure transformation emerges very slowly, so it is possible to observe the various mixed modes representing the superposition of finite-wavelength and longwave modes.The example of such a situation is presented in Figure 6.The intensity of finite-wavelength solutions near the threshold grows according to the square root law, i.e., the excitation of short-wave convection occurs through supercritical bifurcation (Figure 5).
The calculations show that in the longwave mode, where the flow covers not only the fluid layer but also the porous medium, the flow intensity is lower than in the finite-wavelength mode, where the flow is concentrated in the fluid.This is caused by the hydraulic resistance of the porous medium skeleton.Moreover, the intensity of longwave solutions grows very slowly with the increase of the Grashof number.This is why already at small supercriticalities with respect to finite-wavelength disturbances, the intensity of the latter surpasses the intensity of the longwave solution (Figure 5).Apparently, the fast loss of stability of the longwave solution is related to this circumstance.The intensity of finite-wavelength solutions near the threshold grows according to the square root law, i.e., the excitation of short-wave convection occurs through supercritical bifurcation (Figure 5).
The calculations show that in the longwave mode, where the flow covers not only the fluid layer but also the porous medium, the flow intensity is lower than in the finite-wavelength mode, where the flow is concentrated in the fluid.This is caused by the hydraulic resistance of the porous medium skeleton.Moreover, the intensity of longwave solutions grows very slowly with the increase of the Grashof number.This is why already at small supercriticalities with respect to finite-wavelength disturbances, the intensity of the latter surpasses the intensity of the longwave solution (Figure 5).Apparently, the fast loss of stability of the longwave solution is related to this circumstance.Let us discuss the interaction of longwave and finite-wavelength disturbances at small supercriticalities.It turns out that the critical value of the Grashof number for the longwave instability mode is smaller than that for the finite-wavelength mode.Because of this, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size: at the Grashof number values larger than 80 after the transient stage, depending on the initial conditions, a 4, 5, or 6 wave regime is established.
The process of the structure transformation emerges very slowly, so it is possible to observe the various mixed modes representing the superposition of finite-wavelength and longwave modes.The example of such a situation is presented in Figure 6.Let us discuss the interaction of longwave and finite-wavelength disturbances at small supercriticalities.It turns out that the critical value of the Grashof number for the longwave instability mode is smaller than that for the finite-wavelength mode.Because of this, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size: at the Grashof number values larger than 80 after the transient stage, depending on the initial conditions, a 4, 5, or 6 wave regime is established.
The process of the structure transformation emerges very slowly, so it is possible to observe the various mixed modes representing the superposition of finite-wavelength and longwave modes.The example of such a situation is presented in Figure 6.In Figure 7, the dependences of the stationary flow intensity on the parameter C (square root of the Grashof number value) for the disturbances with various values of k are described.On the vertical axis, the squared maximal value of the stream function in the fluid is plotted.
Critical values of the parameter C , obtained by the extrapolation of the amplitude curves on the zero value of the stream function, are in good agreement with the results of linear theory [1].

Conclusions
The calculations show that the finite-wavelength stationary regimes correspond to the convective flow concentrated in the fluid layer and the longwave stationary regimes of the flow occupying both layers.For the parameter values selected for the calculations, the critical value of the Grashof number for the longwave instability mode is smaller than that for the finite-wavelength mode.Because of that, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size.In our opinion, the reason for that is the following.Due to the hydraulic resistance of the porous matrix, the intensity of flow for the In this figure, one of the phases of the long-term transient process of the merging neighboring cells is shown.As one can see, in the central part of a cavity, two vortices of the same sign gradually merge into one vortex (by which point the central vortex has already practically disappeared).
In Figure 7, the dependences of the stationary flow intensity on the parameter C (square root of the Grashof number value) for the disturbances with various values of k are described.On the vertical axis, the squared maximal value of the stream function in the fluid is plotted.In this figure, one of the phases of the long-term transient process of the merging neighboring cells is shown.As one can see, in the central part of a cavity, two vortices of the same sign gradually merge into one vortex (by which point the central vortex has already practically disappeared).
In Figure 7, the dependences of the stationary flow intensity on the parameter C (square root of the Grashof number value) for the disturbances with various values of k are described.On the vertical axis, the squared maximal value of the stream function in the fluid is plotted.
Critical values of the parameter C , obtained by the extrapolation of the amplitude curves on the zero value of the stream function, are in good agreement with the results of linear theory [1].

Conclusions
The calculations show that the finite-wavelength stationary regimes correspond to the convective flow concentrated in the fluid layer and the longwave stationary regimes of the flow occupying both layers.For the parameter values selected for the calculations, the critical value of the Grashof number for the longwave instability mode is smaller than that for the finite-wavelength mode.Because of that, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size.In our opinion, the reason for that is the following.Due to the hydraulic resistance of the porous matrix, the intensity of flow for the Critical values of the parameter C, obtained by the extrapolation of the amplitude curves on the zero value of the stream function, are in good agreement with the results of linear theory [1].

Conclusions
The calculations show that the finite-wavelength stationary regimes correspond to the convective flow concentrated in the fluid layer and the longwave stationary regimes of the flow occupying both layers.For the parameter values selected for the calculations, the critical value of the Grashof number for the longwave instability mode is smaller than that for the finite-wavelength mode.Because of that, at very small supercriticalities the finite-wavelength disturbances fade and the longwave mode is realized.However, as the calculations show, as soon as the Grashof number values become higher than the critical value for finite-wavelength disturbances, the longwave mode loses its stability and the flow is broken into vortices of the small size.In our opinion, the reason for that is the following.Due to the hydraulic resistance of the porous matrix, the intensity of flow for the longwave mode is much lower than that for the finite-wavelength mode.Moreover, it grows very slowly with the increase of the Grashof number value.Because of that, already at the Rayleigh number value, which is just slightly higher than the instability threshold to finite-wavelength disturbances, the flow intensity of the finite-wavelength mode surpasses the intensity of the longwave mode, which results in the fast loss of stability of the longwave regime.

Figure 2 .
Figure 2. Neutral curves for 10 h = and different values of Da .

Figure 2 .
Figure 2. Neutral curves for h = 10 and different values of Da.

Figure 2 .
Figure 2. Neutral curves for 10 h = and different values of Da .

Fluids 2017, 2 Figure 4 .
Figure 4. Streamlines of the longwave regime (only the upper part of the physical domain is shown).

Figure 4 .
Figure 4. Streamlines of the longwave regime (only the upper part of the physical domain is shown).

Fluids 2017, 2 , 39 6 of 8 Figure 4 .
Figure 4. Streamlines of the longwave regime (only the upper part of the physical domain is shown).

Figure 6 .
Figure 6.Mixed regime: the superposition of finite-wavelength and longwave modes.

Figure 7 .
Figure 7.The dependences of stationary flow intensity the parameter C for various values of k :

Figure 7 .
Figure 7.The dependences of stationary flow intensity on the parameter C for various values of k :