Boundary Effect on Marangoni Convection in a Variable Viscosity Fluid Layer

The onset of Marangoni convection in a horizontal fluid layer with a free surface overlying a solid layer heated from below is studied. Problem is focused on the effect of the solid layer depth or its conductivity. The viscosity group, Rv, Biot number, Bi, depth ratio, dr and conductivity ratio, kr, are significant on determining the critical Marangoni number Mc with the corresponding critical wavenumber αc. The characteristics problem is solved numerically. Results show that the temperaturedependent viscosity destabilizes the fluid system but it behaves oppositely when a higher relative thermal conductivity ratio or higher depth ratio is taken.


INTRODUCTION
This study is focused on the instability of Marangoni convection which is induced by a surface tension gradient. The Marangoni instability arises whenever the temperature gradient across the layer exceeds a certain critical value. Realistically, all fluid posses a temperature-dependent viscosity which influences heat transport and the spatical structure of the fluid. The first theoretical study was done by [1] where he suggested there exists a surface tension gradient when he observe a polygonal cellular patterns appear in a paint layers even the paint is on the underside of a plane, horizontal surface. The principle of the exchange of stability has been numerically verified for Pearson's model [2] .
Effect of viscosity plays a major role in this study. In many fluids with large Prandtl number and in particular, in some oils, the fluids posses a temperaturedependent viscosity as viscosity is more sensitive to temperature variations than heat capacity and thermal conductivity and the effects are important on the onset of convection. A variable viscosity was introduced by [3] and starting from his study, effects of variable viscosity is discussed in some problems of the Rayleigh-Bénard [4][5][6][7] , Marangoni [8][9][10][11][12] and also in Bénard-Marangoni instabilities [13][14] . The viscosity variation across the layer can be sufficiently large which are usually well described by exponential laws. Slatchev and Ouzounov [11] studied the steady Marangoni convection under the exponential law of viscosity dependence on temperature with undeformable surface and in microgravity where the Marangoni number increases slightly with viscosity parameter, N from the critical Marangoni number, M c at the free surface and after reaching a maximum at N≈0.7 decreases becoming even smaller than the critical Marangoni number, M c at the free surface. Kalitzora-Kurteva et al. [15] also treated the problem but with deformable surface for conducting and insulating cases. In the case of conducting, the onset of convection exists at the short wavelength (α ≠ 0) and at the long wavelength (α → 0). However, only one mode exists in the insulating case. It shown that the critical values of the Marangoni, M and wavenumbers, α depend mainly on the viscosity variation and weakly to surface deformation.
Later, as density may not be neglected, kinematic viscosity is taken into consideration instead of the dynamic viscosity. Kozoukharova and Roze [12] determined the influence of variable viscosity effect and surface deformation on the convective threshold for the primary steady and oscillatory Bénard-Marangoni in a fluid layer and show that the stability threshold for the short wavelength mode (α ≠ 0) depends strongly on the viscosity variation while the long wavelength (α → 0) is determined by the surface deformation.
However, rapid development of modern techniques in the recent past has posed challenges in studying convective instability problems in much more complicated two and multilayer fluid dynamical systems. It was studied theoretically and experimentally by considering multilayer of fluid, or a fluid layer separated at the middle or bounded from the above or below by a slab.
Scriven and Sternling [16] consider a two layer fluid model in which each layer has infinite depth and examine the local behavior of the system near the interface. They consider the solutocapillary case with mass transfer in either direction; here surface-tension variations are caused by surface-active materials rather than by heat. They allow C r ≠ 0, so that interfacial deformation is permissible and find profound differences in behavior from the Pearson's model. In fact they find that M c = 0, so that the system is always unstable. This zero critical value of M occurs for long wavelength, α → 0. Smith [17] rationalizes this dilemma by reconsidering the one layer model but allowing surface deformation and gravity. He neglects buoyancy, but since the basic state has a hydrostatic pressure field through which the interface deforms, gravity waves are generated. These gravity wave stabilize the long wave instabilities found by [16], and in most practical situations they regain the Pearson [1] result.
Even though single layer systems and double layer systems heated from below have received a great deal of attention in the past, there have been very few studies related to the thermal instability and heat transfer phenomena in a system with more than two layers. Lienhard [18] approach gives exact solutions of the stability problem and the numerical solution presented is much simpler than previous multilayer solutions where he did obtained the critical Rayleigh numbers for the case of three and four fluid layers separated by equally spaced identical midlayers of various thicknesses and conductivities with isothermal outer walls and for the symmetric two-layer problem with outer walls of finite thermal conductivity. Later, Yang [19] study the onset of a system consisting of multilayer fluids separated by solid partitions and the consequent heat transfer increase due to the fluid motion. The solid partitions are found to have optimum efficiency in suppressing fluid motion when they are uniformly distributed. The equivalent Rayleigh number is shown to be the weighted average of the product of the thermal conductivity and the height of each layer.
Yang [13] consider the lower boundary to be a solid plate where it is a perfect insulating boundary condition for thermal disturbances which is difficult compared to conducting boundary condition. It is found that the solid plate with a higher thermal conductivity tends to stabilize the system. The role of the plate thickness is minor in most of the Bénard-Marangoni experiments, while the conductivity of the plate has a significant impact on the stability of the system. Char and Chen [20] focused on Bénard-Marangoni instability with a boundary slab of finite conductivity. They solved the problem numerically and later compared to the asymptotic of the long wavelength. It shown that the critical Rayleigh number, R c increases with thickness of solid layer to the thickness of fluid and thermal conductivity of solid layer to the thermal conductivity of fluid but decrease with Γ = M/R provided the viscosity parameter, B is large. The effect of viscosity variation to the thermal conductivity and thickness of the boundary slab is also discussed in detail. When the viscosity parameter, B is small, it will raise the critical Rayleigh number, R c provided the Biot number, B i is small but behave oppositely when the viscosity parameter, B is large.
In this study, we consider the onset of Marangoni convection in a horizontal fluid layer with the influences of the variable viscosity and the solid plate at the bottom surface. The problem has been solved numerically to obtain a detail description of the marginal stability curves for the onset of Marangoni convection.

FORMULATION OF THE PROBLEM
We examine a horizontal fluid layer of depth d f which is bounded on top by a deformably free surface and at the bottom is subject to a fixed heat flux. The fluid layer is overlying a solid layer of depth d s . We chose rectangular axes with the x and y axes in the plane of the rigid lower boundary and the z axes vertically upwards, so that the bottom surface is given by z = −d s , the solid-fluid interface is given by z = 0 and in the undisturbed state the free surface is located at z = d f . The geometry of the problem are shown in Fig. 1. The density of the fluid is subjected to Newtonian's law: where T is the temperature of the fluid, ρ 0 is its value at a reference temperature T 0 and a is the positive coefficient of the thermal fluid expansion. Following Kozhoukharova and Rozé [14] we consider the kinematic viscosity to be temperature dependent. Here, a linear law for the kinematic viscosity is selected where ν 0 is viscosity at the reference temperature T 0 and 0 T / T | ζ = ∂ν ∂ is assumed constant. The surface tension, σ is assumed to be a linear function of the temperature, 0 0 where σ 0 is the value of σ at the temperature T 0 and the constant γ is positive for most fluids. The fluid is assumed to be an incompressible fluid with variations of viscosity satisfying the continuity equation together with the momentum and the heat equation. These equations are, respectively where i U (i x, y, z) = are the velocity components, T f is the temperature of the fluid, ρ 0 is its density value at a reference temperature T 0 , µ = ρ 0 ν is the dynamic viscosity, κ f is the fluid thermal diffusivity and the pressure inside the fluid is denoted by p. Since our study is focused only for surface tension effect, the buoyancy forces are neglected in Eq. 5.
For the solid layer, the energy equation takes the form where κ s is the solid thermal diffusivity and T s is the temperature of the solid.
When motion occurs, the upper free surface of the layer will be deformable with its position at At the free surface, we have the kinematic condition, together with the conditions of continuity for the normal and tangential stresses, and the heat transfer balance subjected to Newton's law is of the form where P g is the gas pressure, H is the mean curvature of the surface, given by that this choice of scaling was chosen for consistency with the study of Kozhoukharova and Rozé [14] .

LINEARISED PROBLEM
We analyze the linear stability of the basic state by seeking perturbed solutions for any quantity in terms of normal modes in the form where α = (α x 2 +α y 2 ) 1/2 is the total horizontal wave number of disturbance and s is a complex growth rate with a real part representing the growth rate of the instability and an imaginary part representing its frequency. Substituting Eq. 13 into 4-12, we obtain the corresponding linearized equations involving only the Z-dependent parts of the perturbations to take the temperature and the Z-components of the velocity denoted by Θ and W respectively, namely; at Z = 1 together with

SOLUTION OF LINEARIZED PROBLEM
Proceeding in the manner of Hashim and Arifin [9] , we seek solutions in the forms where the exponent ξ and the complex constant A and C are to be determines. Substituting these forms into the Eq. 14 and 15 and eliminating-A and C, we obtain a sixth-order algebraic Eq. For ξ, namely with six different distinct roots which we denote by ξ 1 ,...., ξ 6 . We can use Eq. 18 to eliminate the free surface deflection.
The dispersion relation between M, α, Cr, Bi and Bo is determined by substituting these solutions into the

RESULTS AND DISCUSSION
The onset of Marangoni convection comprising an incompressible fluid with variable viscosity overlying by a solid layer is investigated numerically. As the fluid is subjected to a uniform heat flux below and above (Bi = 0), the critical wave number (α c ) is vanishing and the critical Marangoni number M c attain a constant value 48, which is the exact value known for the case of single fluid layer [1] . In each case investigated in this study, we can identify the critical minima of the marginal stability curves in the (α, M) plane which we denote by M c with corresponding critical wave number α c and hence determine the ranges of M in which all disturbances are stable and those in which unstable disturbances exist. In all cases investigated, we use the value of the viscosity group, R v similar to [14] where the sign of R v depends on the type of luid. If the sign of R v is positive (negative), then the kinematic viscosity is an increasing (decreasing) function of the temperature. For example, the kinematic viscosity of silicon oil decreases when the temperature increase and its give the viscosity group, R v = −0.5 [14] . From a physical perspective, the onset of convection is governed by a sublayer that is more unstable than the full layer if the kinematic viscosity is increase. If R ν = 0, the viscosity of the fluid layer becomes uniform. The critical Marangoni number and wave numbers obtained for different value of R v and d r when k r = 1, C r = 0 for the case Bi = 0 and Bi = 2 are presented in Table 1. From the table we note that increase the value of R v , the critical Marangoni number decreases for both cases, Bi = 0 and Bi = 2.
In general, the critical Marangoni number M c decreases as viscosity group, R v increases and the viscosity group are destabilizing effect to make the system more unstable. For case Bi = 0, with the larger depth ratio, the global minimum occurs at short wavelength (α ≠ 0) and the critical Marangoni number increases. The critical Marangoni number M c increases with d r and the thicker solid layer are clearly a stabilizing factor to make the system more stable because its might store more thermal energy. Figure 2 shows the marginal stability curves of the Marangoni number M with wave number α for different values of d r when k r = 1 and Bi = 0. For the case shown in Fig. 2, the critical Marangoni number, M c increases with the depth ratio, d r and no longer attains at a long wavelength (α → 0). In order to show the contributions of viscosity group, in Fig. 3 and 4, we have reveals that as the viscosity group R v increase, the marginal stability curves shift downwards, hence shows that the viscosity group destabilizes the no-motion state for all wave number for both Bi = 0 and Bi = 2 when k r = 1. As mentioned above, an increase the viscosity group R v , the convective instabilities occurring in the sublayer rather than full layer. Figure 4 shows the marginal stability curves of the Marangoni number, M with wave number, α for a range of values of the viscosity group, R v when k r = 1 and d r = 10 in the case of Bi = 0 and Bi = 2. We can see that the larger ratio of depth, d r = 10 will destabilize the system at a for different values of the thermal conductivity, k r short wavelength (α ≠ 0) where the critical Marangoni number no longer attains at long wavelength (α → 0) for Bi = 0. Figure 5 shows directly the relation between the depth ratio d r and the critical Marangoni number, M c . Two regions are distinguished in this graph. With the depth ratio d r , the critical Marangoni number, M c increases initially, then starts up to a peak value and maintain slightly the same value, depending on the viscosity, R v . When d r → 0, M c becomes lower which conclude that the fluid is more unstable compared to a larger d r . Figure 6-9 indicates the study of thermal conductivity influences to stability curves and is plotted for different viscosity groups. The result is very similar to the earlier case where the viscosity group always has a destabilizing effect on the stability of Marangoni convection. For Bi = 0, the curve have a single global minima at the short wavelength (α → 0). When k r → ∞, M c will increase and prove that an increase of k r will make the fluid become more stable. This is due to the thermal disturbances are easily dissipated deep into the solid layer.

CONCLUSION
Viscosity are clearly a destabilizing factor where we found that an increase of parameter R v causes a decrease of the critical value of the Marangoni number. By increasing the Biot number, the thermal disturbance can easily dissipate into the ambient surrounding and hence lead to a more stable system. When d r → ∞ and kr → ∞, the critical Marangoni number, M c will also increases. In other words, the fluid becomes stable when a smaller R v or larger Bi, d r , k r is taken.