Rayleigh-Bénard-Marangoni Instability in a Self-rewetting Fluid Layer Overlying a Porous Medium

Rayleigh-Bénard-Marangoni instability in a bilayer system of self-rewetting fluid overlying a porous medium is investigated. The upper surface of the fluid layer is assumed non-deformable and a constant temperature gradient is imposed. Unlike in previous works, self-rewetting fluid is considered in this paper, whose surface tension is a quadratic function of temperature. As a result, it would reinforce or suppress the instability caused by buoyancy in different scenarios. Thorough study is carried out and the influence of depth ratio h, Rayleigh number Ra and Biot number Bi on instability of the system is discussed. It was found that the marginal curve is bimodal in certain combination of parameters and mode transition of convection is also found.


Introduction
Marangoni instability, which is originated from a gradient of surface tension, has long been an important topic. For a layer of fluid between horizontal planes heated from below, Bénard (1900) was the first to perform quantitative experiments and observe a pattern of cellular convective motion. It was afterwards analyzed by Lord Rayleigh (1916), attributing such phenomenon to the action of buoyancy forces. Nonetheless, Pearson (1958) later investigated the problem, proposing an alternative explanation for such cellular convective motion that it could also be triggered by surface tension forces. The critical Marangoni number predicted in Pearson's analysis coincides with experiment results. Nield (1964) combined surface tension and buoyancy effects in his study and found out that they reinforce one another and are tightly coupled. Layer of finite extent is consider by Li et al. (2012) and multi-cellular flow structure is perdicted. In recent years, several papers studied the thermocapillary migration of droplets Hu 2012, 2013;Yin et al. 2008).
The study of instability of a bilayer system of fluid and porous layer started by Nield (1977). In this paper, the onset of convection in a fluid layer overlying a layer of porous medium was investigated, in which the upper surface is assumed deformable. Chen (1988, 1989) performed linear stability analysis and experiments on the onset of finger convection in similar system, whose free surface is considered non-deformable. In their study, Darcy's law and Beavers-Joseph condition (Beavers and Joseph 1967) are applied. It was found that instability of the system is bimodal and responds strongly to the variation of depth ratio h. Zhao et al. (2010) investigated thermal effects on Rayleigh-Bénard-Marangoni instability in a similar system, pointing out that a complicated phenomena is overlooked by previous works. Through linear analysis, they pointed out that heat transfer at free surface can lead to mode transition of convection, which only take place in a system with a forceful surface tension effect. More work is done on similar system (Straughan 2001;Samanta 2017).
Different from the fluid in papers mentioned above, nonlinearity in relationship between surface tension and temperature can be found in fluids like aqueous solutions of higher alcohols (Vochten and Petre 1973;ABE 2006), which may show minima at certain temperature. Fluid which exhibits such behavior is often called self-rewetting fluid. Cloot and Lebon (1986) investigated instability of fluid with nonlinear temperature-dependent surface tension and introduce second-order Marangoni number. In their study, error caused by omitting the nonlinear term is also investigated.
To the best of our knowledge, influence of the peculiar characteristic of self-rewetting fluid on the instability of a bilayer system is yet studied. In this paper, we investigated the Rayleigh-Bénard-Marangoni instability in a layer of self-rewetting fluid overlying a porous medium, of which the upper surface is assumed non-deformable. This paper is organized as follows. It begins with mathematical model in "Mathematical Model", in which governing equations and the linearized perturbed system are given. In "Results and Discussion", the influences of depth ratio h, Rayleigh number Ra and Biot number Bi on instability of the system are discussed in detail. A conclusion remark is made in the end.

Mathematical Model
The bilayer system, of which the sketch is shown in Fig. 1, comprises a layer of self-rewetting fluid and a layer of porous medium saturated with the same fluid. The system is considered two-dimensional and infinite in the horizontal direction x. A uniform vertical temperature gradient T∕ z = −b(b > 0) is imposed on the system. The surfacetension exerted at the upper surface is considered to be quadratically dependent on temperature: in which, m represents the minimum surface tension at temperature T min and is a positive parameter.
Employing the Boussinesq approximation of the Navier-Stokes equation with constant dynamic viscosity , kinematic viscosity , thermal diffusivity and coefficient of thermal expansion of the fluid , the governing equations are as follows. At the fluid-porous interface z = 0 , normal velocity, temperature, heat flux, normal momentum are assumed continuous and Beavers-Joseph condition (Beavers and Joseph 1967) is applied.
On the lower surface z = −H m , the boundary is considered to be rigid and perfectly heat conducting. In these equations above, the subscript l, m represent the fluid layer and the porous layer, respectively. g denotes the gravitational acceleration, the porosity of the porous layer, K the permeability, c the specific heat capacity, the thermal conductivity, the convective heat transfer coefficient and the Beavers-Joseph constant. Specially, ( s c s ) represents the heat capacity of the porous medium.
We scale the coordinate, velocity , pressure p, l , respectively. Perturbations of velocities, pressure and temperature are introduced and the equations are linearized thereafter. In order to analyze arbitrary disturbance in terms of normal modes, we suppose that the perturbations w ′ and ′ have the forms where W l , W m , Θ l , Θ m represent the dependence of amplitude on z coordinate, k is the wave-number of the perturbation and is the complex growth rate. Equations 2-7 now become For the fluid layer For the porous layer Boundary conditions are expressed as the following form. On the upper surface z = 1, At the fluid-porous interface z = 0, In the equations above, the dimensionless parameters are listed: where Ra, Pr, , Bi, Ma are Rayleigh number, Prandtl number, Darcy number, Biot number and Marangoni number, respectively. For convenience, we introduce the ratio of heat capacity per unit volume G m = ( c) m ∕ 0 c l , the thermal diffusivity ratio X = l ∕ m and f = T min − T top ∕ΔT . From these equations and boundary conditions, an eigenvalue problem is formulated, which is solved with the Chebyshev numerical method. For consistency with previous

Effect of Surface Tension
The unique characteristic of self-rewetting fluid is the nonlinear dependence of its surface tension to temperature. Unlike fluid in normal Marangoni-Rayleigh-Bénard convection, whose surface tension is supposed to be a monotonically linearly decreasing function of temperature in most cases, self-rewetting fluid could shift from one kind of behavior to another which is absolutely different as temperature varies. In the case when surface tension of self-rewetting fluid decreases as temperature rises, surface tractions tends to encourage the circulation as shown in Fig. 2(a). On the contrary, surface tractions would play a suppressing role when surface tension is a increasing function of temperature as in Fig. 2(b). In Fig. 2, arrows represent the direction along which the surface tension increases, which is also the direction of surface tension forces. In this paper, parameter f in Eq. 26 measures whether the temperature at upper surface is higher, equal or smaller than T m . f > 0 represents the case when temperature is smaller than T min , in which surface tension would decrease as temperature rises. While f > 0 stands for the spurring case, f < 0 would be a signature for the suppressing scenario. Lots of works have been done on the spurring case. Therefore, we focus on the suppressing case and investigate its instability.

Results for Varying Depth Ratio h
The influence of depth ratio h on stability of the system is studied and the results are discussed in this section. The marginal curves are shown in Fig. 3 for various depth ratio h.As surface tension being a suppressing factor, it is worth mentioning that in the area below the curve, the system is unstable and the area above is the contrary case. When depth ration h is smaller than about 4.66, only one mode is observed and the short-wave instability is dominant. As depth ratio h increases, long-wave mode emerges, which is stimulated by perturbation with a comparatively smaller wavenumber. Streamline patterns for both modes are shown in Figs. 4(a), (b), 5(a) and (b), in which the dashed line indicates the fluid-porous interface. For both long-wave mode and short-wave mode, multi-cell structure exists with a rather smaller depth ratio h. As h increases, the convection tends to shift from coupling mode to decoupled mode.
Variation of critical Marangoni number −fMa c is also investigated, as shown in Fig. 6(a). Short-wave instability exists for small depth ratio h and the critical Marangoni number is extremely large, indicating the Marangoni effect would be far from strong enough to stabilize the system and  convection would most likely occur. However, as depth ratio h increases, short-wave mode gets suppressed rapidly. At about h = 4.66 , long-wave mode emerges and yet short-wave mode is still the dominant instability of the system. As depth ratio h continues to increase, short-wave mode gets further suppressed and long-wave mode become dominant at about h = 6.27 . Surface tension soon loses its ability to maintain stability of the system and long-wave instability gets intensified rapidly as h increases. Tough short-wave mode gets even weaker, it continues to exits up to about h = 13.5 . It is noticed that critical wavenumber of short-wave mode is a monotonically linearly decreasing function of depth ratio h and yet critical wavenumber of long-wave mode reaches its minimum at around 6.3.

Results for Varying Rayleigh Number Ra
The effect Rayleigh number Ra has on stability of the system is also investigated. The marginal curves shown in Fig. 7 correspond to cases where Bi = 0 and h = 1 with Raleigh number Ra as variable.
Unlike the previous result on varying depth ratio h, for this combination of parameters, the marginal curve has only one modal. Nonetheless, the system responds to the change of Rayleigh number Ra quickly.
As shown in Fig. 8, the critical Marangoni number −fMa c first increases in an exponential way for Rayleigh number between 2000 and 4000 then surges for larger value of Ra. The variation of maximum and minimum wavenumber of the marginal curve is also examined and certain scaling rules appears to exist for large Rayleigh number Ra, as shown in Fig. 9.

Results for Varying Biot Number Bi
We further investigated the influence of Biot number Bi on stability of the system. Fig. 10(a) and (b) shows marginal curves when h = 2 and h = 5 with Ra = 10000 in both cases.
For insulating case where Bi = 0 , bimodal characteristic is found in h = 5 case and this peculiarity soon vanishes for larger Bi. This shows that certain extent of heat insulation on free surface is essential for long-wave instability. Though marginal curves respond to the increment of Biot number Bi mostly the same way for different depth ratio h, one characteristic is worth noticing. The marginal curves have two intersects. The variation of Biot number Bi has no influence on these points. Despite the fact that they vary as depth ratio h and Rayleigh number Ra changes, the influence is limited.
Although the form of heat transport on free surface is vital to the existence of long-wave mode, its effect on short-wave mode is rather limited. The critical Marangoni number −fMa depends strongly on Biot number Bi. Nonetheless, as the Biot number Bi increases, the critical wavenumber varies little and approaches a value about 3.37 as shown in Fig 11.

Perturbation Energy Analysis
To identify the main instability mechanism, perturbation energy analysis is carried out. Multiplying the linearized disturbance equations by the complex conjugate velocity, we can obtain the perturbation kinetic energy balance equations as follows.
Adding these four equations, integrating across the layer using boundary conditions and taking the real parts gives the following result.
where Each term in Eq. 39 has its physical meaning: E k is the kinetic energy of perturbations, E p is the energy contribution due to compressibility effects on the fluid-porous interface, E vis is the productions as a result of viscous dissipation, E D is the work of Darcy drag and E Ra is the energy contribution due to buoyancy. Specially, E Ma is originated from Marangoni effect on the surface. Since the perturbation kinetic energy E k is positively defined, for the base flow to be unstable, i.e. R > 0 , the sum of terms in the right-hand side of Eq. 39 has the be positive. By investigating the sign and magnitude of each term in the right-hand size, the role each physical effect plays in the destabilization of the base flow and their strength can be identified. It's worth noting that effect of Darcy drag is always a stabilizing factor since E D being negatively defined.
We normalize the perturbation kinetic energy balance equation Eq. 39 so that E k = 1 . Perturbation energy balance of critical states shown in Fig. 6(a) and (b) are investigated, of which the results are shown in Fig. 12(a) and (b). E p has limited contribution to either the stabilization or destabilization of the system. Meanwhile, in both long-wave and shortwave mode, viscous dissipation are the major stabilizing factor in the system. On the other hand, Darcy drag only has a significant effect on stabilizing the system in long-wave mode.
Given that E Ma and E Ra are positive, Marangoni effect and buoyancy are identified as the mechanism which destabilize the system. It is worth pointing out that there is a competitive relationship between these two mechanism. In long-wave mode, buoyancy first contributes the most to the destabilization of the system, then Marangoni effect takes over. In short-wave mode, surface tension is the major contributing mechanism for smaller h and is the weaker one for larger h.

Conclusion
In this paper, we investigated the Marangoni-Rayleigh-Bénard instability in a self-rewetting fluid layer overlying a porous medium subjected to a constant temperature gradient. In contrast to most previous papers where surface tension is assumed to be a linearly decreasing function of temperature, a quadratic relationship is considered. Modes of convection in the system is studied with linear stability analysis. Different effects of surface tension is analyzed showing that not only can surface traction encourage the spur of convection, it can also suppress instability in the system. Effect of depth ratio h, Rayleigh number Ra, Biot number Bi on stability of the system are examined in detail. Bimodal characteristic is observed in some combination of parameters and convection can only be found in fluid layer. Though surface tension is able to play a suppressing role on instability of the system, its effect is rather limited. While short-wave mode dominates when depth ratio h is small, it gets suppressed for larger h. Meanwhile, long-wave mode emerges and becomes the dominant instability for h > 6.27 . Though the bimodal characteristic responds sensitively to h, Biot number Bi is also essential for its existence. However, Bi has limited effect on short-wave mode. Rayleigh number Ra does not contribute to the emergence of longwave mode, yet having its significance on short-wave instability. Through perturbation energy analysis, Marangoni effect and buoyancy are identified as the destabilizing mechanism in the system.