Instabilities of Thin Viscous Liquid Film Flowing down a Uniformly Heated Inclined Plane

Convective instabilities driven by buoancy (Rayleigh-Be ́nard) and /or thermocapillary (Marangoni) effects are very important and fascinating problems in fluid mechanics. In thin liquid film buoyancy stabilizing mechanisms are negligible while thermocapillary instabilities become dominant.

In industry the Marangoni effects have been utilized in a process known as Marangoni drying    [1, 2]. The main advantage of this drying process is the processing speed. It has an extensive use in semiconductor industry to dry silicon wafers. An extensive review of literature on non-isothermal free surface flows, in various context, can be found in Oron et al. [3], Nepomnyashchy et al. [4], Colinet et al. [5] and Velarde and Zeytounian [6].

In the context of thin film, Goussis and Kelly [7] first investigated the thermocapillary instability in a liquid film falling down on an inclined uniformly heated plane. They performed linear stability analysis based on Orr-Sommerfeld and linearized energy equations. They examine the occurrence of thermocapillarity, surface wave instability and detect three different mechanismsby which energy is transferred to the disturbance. Depending on any of these three mechanisms is the dominant one; the instability will first assume the form of either transverse waves or longitudinal rolls. Joo et al. [8] investigated the long wave instability of a uniform film of volatile viscous liquid falling down a uniformly heated inclined plane for two-dimensional disturbances taking into account the surface tension, thermocapillary and evaporation effects. That study is actually an effective extension of Benney’s [9] long wave expansion for non-isothermal film flow. However, in the relevant evolution equation, known as ‘Benney equation with Marangoni effect’ in the literature, the terms depending on the Marangoni number are of second order in comparison with those due to gravity and mean flow. They applied the linear and nonlinear stability theory to analyze the competition among different types of instabilities and obtained similar results as those of Goussis and Kelly [7].

A considerable amount of  study has been focused on the long wave instability for a liquid film flowing down a uniformly heated plane (see for example Boos and Thess [10]; Kalliadasis et al.[11]; Trevelyan et al. [12]; Scheid et al. [13]; Ruyer-Quil et al. [14]; Scheid et al. [15]) and also for no uniformly heated plane (see for example Miladinova et al. [16, 17]; Mukhopadhyay and Mukhopadhyay [18]) and related studies that extend these analyses (see for example Yeo et al. [19]; Saprykin et al. [20]). Out of these studies Scheid et al. [13] examined the validity domain of the Benney equation with Marangoni effects. They find that the thermocapillary can strongly reduce the validity domain of the Benney equation. They also ensure that the Benney equation still remains helpful to study thin film flows with various physical effects as it allows describing the film evolution by a single equation with various relevant parameters, in spite of some limitations in its validity domain.

In the present study we have derived a Benney type equation with Marangoni effect and specially discussed the effect of Biot number which describes heat transfer at the free surfaceon the thin film flow.

The rest of this study is organized as follows: Section 2 describes the formulation of the problem and derivation of the free surface equation. Section 3 is devoted to linear and weakly nonlinear stability analysis of the problem considered followed by concluding remarks.

2.Formulation of the Problem

  A thin viscous liquid film of density , dynamic viscosity μ, thermal conductivity and surface tension  flowing down a uniformly heated plane of constant temperature distribution is considered. The thin liquid film is bounded by ambient air of temperature  and pressure  and its inclination to the horizon is . The dynamic influence of the air above the film is ignored. Also, the plane is assumed to be a perfect heat conductor. The temperature difference between the plane and the bounded air induces a thermocapillary Marangoni effect which affects the free surface and the flow of the liquid. The surface tension of the liquid is assumed to vary linearly, within the range of small temperature difference as:

                                                  (1)

Where is the surface tension at , which is taken as the reference temperature and , the thermal surface tension coefficient is a positive constant for most common fluids. The assumption of linear variation of surface tension with temperature is much compatible with the experimental data for many fluids (O’ Brien [2]). The other properties of the liquid such as density , kinematic viscosity and thermal diffusivity  with  the specific heat at constant pressure etc. are assumed to be constant. Although it is well-known that the temperature variation in the liquid can cause dramatic changes in the above mentioned properties, so approximation can be made depending on the type of the problem being examined. The heat transfer coefficient from the liquid to the air is denoted by . A rectangular coordinate system is chosen where the axis is along the inclined plane and axis is normal to the inclined plane. The main objective of this study is to examine the effect of heat transfer at the free surface on the stability characteristics of the liquid film.

The governing equations are the conservation of mass, momentum for the fluid flow and energy equation for the temperature field

Where  is the velocity vector, g=  is the gravitational acceleration and T denotes the absolute temperature.

The boundary conditions are the usual no – slip condition and the constant temperature distribution as the thermal boundary condition on the inclined plane  :

V=0,                                                                                                                 

At the free surface , dynamic and kinematic conditions along with Newton’s law of cooling as the thermal boundary condition are,

Where n and t are unit vectors along the normal (outward pointing) and tangential to the interface respectively and  is the stress tensor with the rate of strain tensor

For writing the problem in dimensionless form, the dimensionless quantities are defined as:

                                                                                                                    (4)

Where we assume l₀ as the characteristic longitudinal length scale whose order may be considered the same as the wave length , the mean film thickness  as the length scale in transverse direction and the Nusselt velocity as the characteristic velocity.

Using the dimensionless quantities (4) in the governing equations and boundary conditions (2-a) - (3-f) we arrive after dropping the asterisk as:

                                                                                                        (5)

                                                                                   (6)

                                                                                                                 (7)

                 (8)

(II) Boundary conditions at the inclined plane (z=0):

                                                                     (9)

(III) Boundary conditions at the free surface

                                                                       (10)

                                  (11)

                                                                                                (12)

                                   (13)

 Where   is the Reynolds number and We  is the weber number; Fr  is the Froude number (which is related with Reynolds number as  is the Prandtl number  is the Marangoni number, Sr  is the surface-tension number, Bi  is the Biot number and  is the aspect ratio.

It is assumed that the Reynolds number is O(1), Weber number is , Marangoni number is of order , Surface-tension number is , the Prandtl number is , and the Biot number is : An estimation of the dimensionless parameters is given in appendix 1.

The complexity of the free surface problem posed by equations (5) -(13) can be overcome by invoking long-wave (lubrication) approximation for . This provides us an opportunity for asymptotic reduction of the governing equations and boundary conditions to a single non-linear partial differential equation of Benney [9] type, in terms of the local film thickness h(x, t). A comprehensive review on long-wave approximation method can be found in Chang [21] and in many other studies. Introducing the stream function , defined by

,                                                                                                  (14)

in the governing equations and the pertinent boundary conditions (5)- (13) and considering the above mentioned orders of the non-dimensional parameters, we obtain

,                                                                                                                                                                (15)

,                               (16)

                           (17)

                                          (18)

                                                                          (19)

                                                                                             (20)

                                      (21)

                                           (22)

Where Pe = RePr is the Peclet number for heat transport that measures the relative importance of conduction over convection, which is of   according to our assumption.

We are now interested in yielding a non-linear evolution equation in terms of non-dimensional film thickness h(x, t), depending on the dimensionless spatial and temporal variables x and t.The equation would represent the spatiotemporal dynamics of the two dimensional liquid filmof a mean thickness h₀. Since the long -wave length modes are the most unstable ones for film flow, expanding the physical quantities , p and T as a power series,

and substituting the above into the governing equations (15)-(17) and the boundary conditions (18)-(21) and then collecting the coefficients of like powers of  the zeorth and the first-order equations are obtained as :

I) Zeorth order equations:

          (23)

                     (24)

                                                                      (25)

II) First-order equations:

                                                                           (26)

                                                                                             (27)

                           (28)

                                                       (29)

                        (30)

                                                                                                                            (31)

                                                               (32)

Now solving the equations (23)-(32) the zeroth and the first-order solutions are obtained as:

(I) Zeorth order solutions:

                           (33)

                                                                                (34)        

           (35)

Let us note that the thermocapillarity influences the longitudinal velocity even at the zeroth order approximation.

(II) First-order solutions

                          (36)

                                               (37)

                                                                               (38)

Where  are given in the appendix 2. Substituting   and     into the kinematic boundary condition (22) we obtain a nonlinear generalized kinematic equation in terms of the surface deformation h (x,t) as :

                                                (39)

and the suffixes denote the differentiation with respect to the corresponding variables,

The equation (39) is Benney equation with Marangoni effect and is similar to the equation that is derived by Joo et al. [8] for uniformly heated plane(in the absence of evaporation effects /intermolecular forces) but it contains some additional terms with coefficients C(h), D(h), E(h), F(h) and G(h) which are responsible for thermocapillary effects only. The other terms are also modified by the thermocapillarity.

Let us point out here that the Biot and Marangoni numbers appear in the equation (39) in such a way that the limit of Mr → 0 or Bi → 0 leads to the evolution equation corresponding to that of isothermal case. A close scrutiny shows that if any one of them tends to be zero, other will also tends to be zero, in this case. The reason is that, for Mr → 0 implies T_w = T_g and so there is no temperature gradient across the interface and consequently no heat transfer  takes place through the interface and so k_g→ 0 implies that Pi → 0. On the other hand, using equations (34) and (37), Bi →0 implies that  which implies that the temperature is everywhere uniform and equal to the wall temperature. Therefore, T_g = T_w, which leads to Mr → 0. This is not true in general because the two parameters represent the different and unrelated physical processes. The Marangoni number may not be zero, even if there exists a temperature gradient across the film, even for vanishing Biot number, as discussed by Mukhopadhyay and Mukhopadhyay [18].

3. Stability Analysis

The variation of the film thickness in the basic flow is found to be very small, so it is reasonable to assume that the local dimensionless film thickness is equal to one. Thus, the dimensionless film thickness in the perturbed state may be written as

                                                                      (40)

Where  is the dimensionless perturbation of the film thickness. Setting the transformation as

                                                                           (41)

and using equations (40) and (41) in equation (39), and retaining the terms up to the second-order fluctuations after dropping the tilde sign can be written as

Where A, B, C, D, E, F, G and their corresponding derivatives, which are denoted by the primes, are evaluated at h = 1. Equation (42) characterizes the behaviour of finite amplitude perturbation on the thin viscous film and it is used to signify the time sapient behaviour of an initially sinusoidal perturbation on the viscous film.

Let us note that the assumption of constant film thickness with long wave perturbation for the basic flow as used in equation (40) is judicious one, for otherwise if we take h(x, t) = h_B(x)+ (x, t) instead of equation (40), where h_B(x)is the time-averaged part of h(x, t), then to obtain equation (42) we have (i)to neglect the derivatives of the averaged film thickness h_B(x)with respect to x in the unsteady equation and thence (ii) to allow h_B(x)→ 1 in the unsteady equation. The first assumption corresponds to the parallel flow approximation and the second leads to the restriction that the unsteady equation (42) is only locally valid.

3.1.Linear stability analysis

Following standard linear stability analysis [18], we are interested to study the linear response for a sinusoidal perturbation of the film by assuming the perturbation of the form

                             (43)

Fig.1 Variation of  for different  and for .

  =5.5in

Where  is the amplitude of the disturbance and c.c represents the complex conjugate. Here the wave number k is real and is the complex frequency. Using equation (43) in linearized part of equation (42), we get the dispersion relation as

                                                                                                                                                                                          (44)

Equating the real and imaginary parts of equation (44), we get

                                                                                                              (45-a)

and

                                                                                                            (45-b)

The onset of instability can be obtained by equating  to zero, which gives the critical condition.

           (46)

The expression (46) has also the same functional form as the one derived by Kalliadasis et al.[11] except the possible difference, as that is calculated by momentum integral method.

The above relation (46) defines the critical Reynolds number

                                                     (47)

Above which the flow loses stability.

Fig.1 depicts the variation of the critical Reynolds number with the Biot number for different Marangoni numbers. The figure reveals that increasing Marangoni number Mr decreases critical Reynolds number i.e. Marangoni number enhances the destabilizing mechanism which is compatible with the results found by Mukhopadhyay and Mukhopadhyay [18]. Also, the figure shows that  decreases with increasing Bi and reaches its minimum value at Bi =1 and then increases with Bi. Thus within the range 0 < Bi < 1 instability increases with Bi. This is due to the fact that in this range of Bi heat conduction in thefluid region is more comparable to that of the surrounding air, as a consequence interface will be excited more and more with increasing the Bi. At Bi = 1,  is minimum or in other words destabilizing effect is maximum. This is due to the fact that for Bi = 1 free surface behaves as a thermally insulated one and therefore no heat transfer takes place through the interface or in other words thermal energy is conserved within the fluid layer and so thermal excitement reaches its optimum stage or we can say that for Bi = 1, there is maximum room for interfacial perturbation.

The stabilizing mechanism of Bi in the range Bi > 1 can be fairly understood as in this case the surrounding air is more conducting than that of the fluid, as a consequence the thermal energy is wiped out quickly from the interface.

Let us draw our attention to the fact when Mr → 0 or Bi → 0,  = (5/6) cot , which is the critical Re, for isothermal viscous film flowing down an inclined plane under the action of gravity only.

The neutral stability condition  gives two relations

                                                                          (48-a)

and

               (48-b)

Fig. 2  Variation of  for different  and for , ,  and .

  =5.5in

Fig. 3  Variation of  for different  and for , , ,  and .

=5.5in

Which correspond to two branches of the neutral curves and the flow instability takes place between them. Further the neutral curves intersect at the bifurcation point Re = k = 0. The wave number of the waves with a maximum growth is obtained from the relationd /dk = 0, which gives .

Fig.2 shows the variation of the neutral curve k_n with Bi for different Marangoni numbers.It is clear from the figure that as Mr increases k_n also increases. Moreover, as Bi increases k_n also increases and reaches its maximum at Bi =1 and then decreases  with increasing the Bi. This result confirms the behaviour of Bi as discussed earlier to explain the Fig.1. Fig.3 depicts that the growth is maximum at Bi =1, as expected.

In the next section we are interested to study the growth of weakly nonlinear waves and therefore we need to consider two or more characteristic times scales for different orders of magnitudes. As a rule, problems of this kind cannot be solved by means of classical perturbation method and the method of matched multiple scale expansions is also insufficient.

3.2. Weakly nonlinear stability analysis

Now we are interested to study those small amplitude waves which develop immediately just after the breakdown of the flat film solution

From the linear stability analysis we see that in the neutral state all modes are neither stable nor unstable and therefore , which corresponds to two branches of the neutral curve and the instability takes place between them. Furthermore, the neutral curve intersect at the bifurcation point . Thus, one  expects in the vicinity of the upper branch of the neutral curve, a thin band of width 1 (say) of unstable modes around a central one  which appears over a time  such as . Therefore, the main features of the system behaviour can be obtained from an analysis of these modes and so the nonlinear evolution of the waves is investigated in the region where .

The above expectation, together with the anticipation that the wave packet may travel at a group velocity of order one, suggest a scale

                                                             (49)

for the multiple scale analysis as outlined in Debnath [22]. Here t and x are fast scales, whereas x₁, t₁ and the rest are slow scales. It is assumed that these variables are mutually independent, so the temporal and spatial derivatives become

                                                                                                                                                                     (50)

Now, the surface  deformation is expressed as

                                                                           (51)

Using equations (49), (50), (51), in equation (42) we get

                                                                                     (52)

Where , ,  etc. are the operators and N₂, N₃ are the nonlinear terms of equation  given in the appendix 3.

In the lowest order of , we have

                                                                                                              (53)

Which has a solution of the form

                          (54)

Where  and c.c. denotes the complex conjugates. It is to be noted here that the above solution given in equation (54) is already obtained in connection with the linear stability analysis except is replaced by , since in the vicinity of the neutral curve , so that the function exp( ) is slowly varying and may be absorbed in In the second order, the perturbation system yields

                                               (55)

Fig. 4  Variation of  for different  and for , , ,  and .

=5.5in

Invoking equation (54)  in equation (55) we have

                                                                          (56)

Where c.c. denotes the complex conjugate and the notations H₂, Q₁ are given by

The function has the functional form , since we consider modal interaction near the neutral curve such that . It is to be noted here that the secular condition leads to the requirement that the wave envelope must propagate at the speed H_2r which turns out to be the group velocity, since . Eliminating secular terms the solution for  is obtained  from equation (56) as

        (57)

Where

Substituting  and in to the equation (52) of third order in  after elimination of its secular term gives

                                               (58)

Where

,         and

Using the transformation

and 

in the equation (58)we get,

  (59)

Which is standard Complex Ginzburg-Landau Equation (CGLE). In a different context CGLE is derived by Dandapat and Samanta [23] from second-order Benney equation.

To find the linear stability of the amplitude equation (58) we first linearize the equation at about and then following the normal mode analysis substitute with complex growth rate  and wave number k, both dimensionless, in the linearized equation, then we get the dispersive relation

                                                    (60)

The disturbance grows with time about  if

                                                                  (61)

To satisfy the above inequality the wave number k must be of the following forms

                                                                                   (62)

for  or

for

In the nonlinear regime the stability characteristics of the perturbation dynamics solely depend on the sign of J₂. If the sign is J₂> 0, then the disturbance amplitude decays withtime in the unstable region. This is the case of a supercritical bifurcation. On the other hand, if the sign is J₂< 0,then the disturbance amplitude grows with time in the stable region.This is the case of a subcritical bifurcation. Careful observation of the fig. 4 shows that  increasing the Marangoni number Mr diminishes the supercritical bifurcation zone and its influence is very strong and significant. It is interesting to note that there is a critical value of theReynolds number Re before which J₂ increases with  increasing the Mr while it behaves just opposite after this critical value of Re.

 4. Conclusion

The stability of a thin viscous film flowing down a uniformly heated plane is investigated by long wave expansion method. Special attention is given to scrutinize the effect of Biot number on the instability mechanism. Based on the above analysis given in sections 3.1 and 3.2 the following conclusion can be  reached.

(1) The linear stability analysis of the evolution equation shows that the Biot number playsa double role; for Bi < 1 it gives destabilizing effect but for Bi > 1 it produces stabilization. At Bi = 1 the instability is maximum. The linear stability analysis also shows that  increasing the  Marangoni number Mr, decreases the critical Reynolds number  i.e. Marangoni number enhances the destabilizing mechanism.

(2) The weakly nonlinear study reveals that the impact of Marangoni number Mr is very strong on bifurcation scenario even for a slight variation.

 Acknowledgements

The authors express their gratitude to the referees for their valuable comments to improve the text. Anandamoy Mukhopadhyay expresses his sincere thanks to The University of Grants Commission, India, for providing necessary grants vide Minor research project No. PSW- 025/11-12 (ERO) dated 03/08/2011 and Sanghasri Mukhopadhyay expresses her sincere thanks to The Department of Science and Technology, Government of India for awarding INSPIRE scholarship, Registration number: [1964/2011], to pursue the higher education and research.

Appendix 1

To make the assumption of the order of the dimensionless parameters realistic, we like to estimate the approximate ranges of the parameters by considering some real physical data. We have calculated the mean film thickness of the film from the relation by setting a fixed Reynolds number Re and . The physical data (collected from Yeo et al. [19]) and the estimated value of the relevant parameters are given in TABLE I andII below.

Convincing from the estimated values of the parameters given in Table II we assume that the Reynolds number is O(1), Weber number is , Marangoni number is of order , Surface-tension number is , the Prandtl number is O(1) and the Biot number is O(1).

Table.I Typical values for the relevant physical constants for water/silicon oil 50 cS

Table.II Typical values for the relevant dimensionless parameters for water/silicon oil 50 cS

Appendix 2

Appendix 3

References

[1] J. Marra, J. A. M. Huethorst, Physical principles of Marangoni drying, Langmuir, 7, 2748-2755 (1991).

[2]           S. B. G. M. O’ Brien, On Marangoni drying : nonlinear kinematic waves in a thin film, J.Fluid Mech., 254, 649-670 (1993).

[3]           A. Oron, Nonlinear dynamics of thin evaporating liquid films subject to internal heat generation, In Fluid Dynamics at Interfaces (ed. W. shyy and R. Narayanan), Cambridge University Press, (1999).

[4]           A. A. Nepomnyashchy, M. G. Velarde, P. Colinet, Interfacial phenomena and convection, Chapman and Hall, (2002).

[5]           P. Colinet, J. C. Legros, M. G. Velarde, Nonlinear dynamics of surface-tension driven instabilities, Wiley VCH, (2001).

[6]           M. G. Velarde, R. Kh. Zeytounian, Interfacial Phenomena and the Marangoni effect, Springer, (2002).

[7]           D. A. Goussis, R. E. Kelly, Surface wave and thermocapillary instabilities in a liquid film flow, J. Fluid Mech., 223, 25-45, (1991) and corrigendum J. Fluid. Mech. 226, 663- (1991).

[8] S. W. Joo, S. H. Davis, S. G. Bankoff, Long-wave instabilities of heated falling films: two-dimensional theory of uniform layers, J. Fluid Mech., 230, 117-146 (1991).

[9]           D. J. Benny, Long waves on liquid films, J. Math. Phys., 45, 150-155 (1966).

[10]         W. Boos, A. Thess, Cascade of structures in long-wavelength Marangoni instability, Phys.Fluids, [1], 1484-1494 (1999).

[11]         S. Kalliadasis, E. A. Demekhin, C. Ruyer-Quil, M. G. Velarde, Thermocapillary instability and wave formation on a film falling down a uniformly heated plane, J. Fluid Mech., 492, 303-338 (2003).

[12]         P. M. J. Trevelyan, S. Kalliadasis, Wave dynamics on a thin-liquid film falling down a heated wall, J. Engineering Mathematics, 50, 177-208 (2004).

[13]         B. Scheid, C. Ruyer-Quil, U. Thiele, O. A. Kabov, J. Legros, P. Colinet, Validity domain of the Benney equation including Marangoni effect for closed and open flows, J. Fluid Mech., 527, 303-335 (2005).

[14]         C. Ruyer-Quil, B. Scheid, S. Kalliadasis, M. G. Velarde, R. Kh. Zeytounian, Thermocapil-lary long waves in a liquid film flow, Part 1 Low-dimensional formulation, J. Fluid Mech., 538, 199-222 (2005).

[15]         B. Scheid, C. Ruyer-Quil, S. Kalliadasis, M. G. Velarde, R. Kh. Zeytounian, Thermocapil-lary long waves in a liquid film flow, Part 2 Linear stability and nonlinear waves, J. Fluid Mech., 538, 223-244 (2005).

[16]         S. Miladinova, S. Slavtchev, G. Lebon, J. C. Legros, Long-wave instabilities of non-uniform heated falling films, J. Fluid Mech., 453, 153-175 (2002).

[17]         S. Miladinova, D. Staykova, G. Lebon, B. Scheid, Effect of non-uniform wall heating on the three-dimensional secondary instability of falling films, Acta Mechanica, 30, 1-13 (2002).

[18]         A. Mukhopadhyay, A. Mukhopadhyay, Nonlinear stability of viscous film flowing down an inclined plane with linear temperature variation, J. Phys. D: Appl. Phys., 40, 1-8 (2007).

[19]         Y. L. Yeo, V. R. Craster, O. K. Matar, Marangoni instability of a thin liquid film resting on a locally heated horizontal wall, Physical Review E, 67, 056315 1-14 (2003).

[20]         S. Saprykin, P. M. J. Trevelyan, R. J. Koopmans, S. Kalliadasis, Free surface thin film flows over uniformly heated topography, Physical Review E, 75, 026306 1-17 (2007).

[21]         H. C. Chang, Wave evolution on a falling film, Annu. Rev. Fluid Mech., 26, 103-136 (1994).

[22] .       L. Debnath, Nonlinear Partial Differential Equations for Scientists and Engineers, Springer, (1997).

[23]         B. S. Dandapat, A. Samanta, Bifurcation analysis of first and second order Benney equa-tions for viscoelastic fluid flowing down a vertical plane, J. Phys. D: Appl. Phys., 41, 095501(9PP) (2008).