Series Solutions for Marangoni Convection on a Vertical Surface

Abstract

The problem of steady, laminar, thermal Marangoni convection flow of Newtonian fluid over a flat surface is investigated. The boundary layer equations for the momentum and energy equations are transformed with the similarity solutions to ODEs to obtain the analytical approximate solutions. The analysis assumes that the temperature variation is a power law function of the location. The approximate solutions to the similarity equations are obtained by exponential series. The effects of the power law exponent and Prandtl number on the velocity and temperature profiles are presented.

1. Introduction

For many metallurgical and materials processing applications in space, it has been evidenced that if a free liquid surface (or liquid-liquid interface) is present, a surface tension gradient of this free surface will cause liquid movement at this free surface which alternatively drives fluid movements in the corresponding phase. This phenomenon is referred to as capillary motion. When a free liquid surface is present, the surface tension variation resulting from the temperature or concentration gradient along the surface can also induce motion in the fluid called solutal capillary, and thermocapillary motion, respectively.

The study of liquid movement resulting from thermocapillarity (or so-called Marangoni) convection is very important for a liquid system either in microgravity or in normal gravity [1]. Under normal gravity, liquid movement is mainly driven by buoyancy force because of the temperature-dependent density while the liquid is exposed to a temperature gradient field. As the size of the liquid system decreases especially having the size decrease in the direction of gravity, the buoyancy effect begins to diminish and the Marangoni effect will then dominate the system as the main driving force for liquid interface movement. In the absence of gravity, Marangoni convection always plays a main role in the determination of the fluid movement because of varying liquid surface tension in a temperature gradient field regardless of the system’s size. It has significance in the processing of materials, especially in small-scale and low-gravity hydrodynamics [2] to investigate the interaction between the heat and mass transfer.

Marangoni convection appears in many industrial processes and space technologies, for example, in the study of semiconductor crystallization processes, in crystal growth melts where the flow produces undesirable effects (see [3–6]) and it occurs around vapor bubbles during nucleation [5]. In several papers authors investigate Marangoni driven boundary layer flow in nanofluids. These fluids can tremendously enhance the heat transfer characteristics of the base fluid and have many industrial applications such as lubricants, heat exchangers, and coolants. Nanofluids are studied when different types of nanoparticles such as oxide ceramics, nitride ceramics, and aluminium, are used [7–9]. Marangoni flow has also significance in welding, semiconductor processing, and other fields of space science. Its mathematical model is studied in [7, 10, 11]. Marangoni boundary layers may occur along the liquid-liquid or liquid-gas interferences.

These phenomena have been investigated for various substances in geometries with flat surfaces by similarity analysis (see [3, 7, 12–16]). Arafune and Hirata [3] presented a similarity analysis for just the velocity profile for Marangoni flow when the surface tension variation is linearly related to the surface position. Christopher and Wang [4] studied Prandtl number effects for Marangoni convection over flat surface and presented approximate analytical solutions for the temperature profile. They showed that the calculated temperature distribution in vapor bubble attached to a surface and in the liquid surrounding the bubble was primarily due to the heat transfer through the vapor rather than in liquid region and the temperature variation along the surface was not linear but could be described by a power-law function [12]. Using the similarity transformation, the governing system of nonlinear partial differential equations are transformed into a pair of similarity nonlinear ordinary differential equations, one for the stream function and one for the temperature. The velocity and temperature distributions can be given by numerically using the Runge-Kutta method [5, 7, 12, 17], analytical approximate solutions can be determined for these problems by using Adomian decomposition method and Padé technique [15, 16, 18, 19] or by power series method [20].

This paper investigates a similarity analysis for Marangoni convection inducing flow over a flat surface due to an imposed temperature gradient. The analysis assumes that the temperature variation is a power law function of the location and the surface tension is assumed to depend on the temperature linearly.

We first present the derivation of the equations and show how the boundary layer approximation leads to the two points boundary value problem and the similarity solutions. The new model, written in terms of stream function and temperature, consists in two strongly coupled ordinary differential equations. Its analytical approximate solutions are represented in terms of exponential series. The influence of various physical parameters on the flow and heat transfer characteristics are discussed.

2. Boundary Layer Equations

Consider the steady laminar boundary layer flow of a viscous Newtonian fluid over a flat surface in the presence of surface tension due to temperature gradient at the wall. Assuming that the surface is impermeable, the surface tension varies linearly with temperature and the interface temperature is a power-law function of the distance along the surface. The governing equations for two-dimensional Navier-Stokes and energy equations describing thermocapillary flows in a liquid layer of infinite extent are considered. The layer is bounded by a horizontal rigid plate from one side and opened from the other one. The rigid boundary is considered as thermally insulated. The physical properties of the liquid are assumed to be constant except the surface tension. This investigation are based on the balance laws of mass, momentum and energy. In the domain 𝑥>0,𝑦>0, these can be written in the form [2]:𝜕𝑢𝑢𝜕𝑢𝜕𝑥+𝜕𝜕𝑣𝑦=0,(2.1)𝜕𝑥+𝑣𝜕𝜕𝑦𝑢=𝐾𝜕²𝑢𝜕𝑦²𝑢𝜕𝑇,(2.2)𝜕𝑥+𝑣𝜕𝜕𝑇𝑦=𝛼^𝑓𝜕²𝑇𝜕𝑦²,(2.3) where the 𝑥 and 𝑦 axes are taken parallel and perpendicular to the surface, 𝑢 and 𝑣 are the parallel and normal velocity components to the surface, respectively, and 𝛼_𝑓 denotes the thermal diffusivity, 𝐾=𝜇/𝜌 is the kinematic viscosity and density ratio of the ambient fluid.

Marangoni effect is incorporated as a boundary condition relating the temperature field to the velocity. The boundary conditions at the surface (at 𝑦=0) are𝜇𝜕𝑢𝜕𝑦||||_𝑦=0=−𝜎_𝑇𝜕𝑇𝜕𝑥|||_𝑦=0𝑇(𝑥,0)𝑣=(𝑥𝑇,(00),=0)0+,𝐴𝑥,_𝑚+1𝑢𝜕(𝑇𝑥,∞)=0,,𝜕𝑦||||_𝑦=∞=0,(2.4) as 𝑦→∞, where 𝜎_𝑇=𝑑𝜎/𝑑𝑇, 𝐴 denotes the temperature gradient coefficient, 𝑚 is a parameter relating to the power law exponent. The case 𝑚=0 refers to a linear profile, 𝑚=1 to the quadratic one. The minimum value of 𝑚 is −1 which corresponds to no temperature variation on the surface and no Marangoni induced flow.

Introducing the stream function 𝜓 by𝑢=𝜕𝜓𝜕𝑦,𝑣=−𝜕𝜕𝜓𝑥,(2.5) equation (2.2) is reduced to𝜕𝜓𝜕𝑦𝜕²𝜓𝜕𝑦𝜕𝑥−𝜕𝜕𝜓𝑥𝜕²𝜓𝜕𝑦²=𝐾𝜕³𝜓𝜕𝑦³.(2.6)

Applying similarity functions𝜓=𝐶₁𝑥^𝑎𝑓(𝜂),Θ=𝑇−𝑇(0,0)𝐴𝑥^𝑚+1,(2.7) and similarity variable 𝜂=𝐶₂𝑥^𝑏𝑦 with 𝐶₁=√³(𝑚+1)𝜇𝜎_𝑇𝐴/𝜌²,𝐶₂=√³(𝑚+1)𝜌𝜎_𝑇𝐴/𝜇²,𝑎=(𝑚+2)/3 and 𝑏=(𝑚−1)/3 one can obtain from the system of partial differential equations (2.5), (2.6) one single ordinary differential equation of the third order𝑓^−2𝑚+13𝑓^2+𝑚+23𝑓𝑓^=0,(2.8) and boundary conditions (2.4) become𝑓(0)=0,𝑓^(0)=−1,𝑓^(∞)=0.(2.9) For (2.3) by the similarity temperature function Θ with the corresponding boundary conditions we get(𝑚+1)𝑓^Θ−𝑚+23𝑓Θ^=1PrΘ^′Θ(0)=1,Θ_,(2.10)(∞)=0,(2.11) where Pr=𝜇/(𝜌𝛼_𝑓) is the Prandtl number. For the dimensionless stream function 𝑓(𝜂) and the temperature field Θ(𝜂), the system (2.8), (2.10) is derived and the primes denote the differentiation with respect to 𝜂.

Now, the velocity components can be expressed by similarity function 𝑓 as follows𝑢(𝑥,𝑦)=𝜕𝜓𝜕𝑦=𝜅²³𝜌𝜇𝑥^(2𝑚+1)/3𝑣(𝑥,𝑦)=−𝜕𝜓𝑓^′(𝜂),𝜕𝑥=−𝜅3𝑥^(𝑚−1)/3(𝑚+2)³𝜇𝜌𝑓(𝜂)+(𝑚−1)𝜅³𝜌𝜇𝑥^(𝑚−1)/3𝑦𝑓^(𝜂),(2.12) where 𝜅=√³(𝑚+1)𝜎_𝑇𝐴/𝜌, 𝜎_𝑇= constant.

It should be noted that 𝑢 and 𝑣 are proportional to 𝑥^(2𝑚+1)/3 and 𝑥^(𝑚−1)/3, respectively. It means, that for 𝑚=−1/2 the velocity component 𝑢 is a constant on the upper surface of the boundary layer. If 𝑚=1 then 𝜂=√³(𝑚+1)𝜌𝜎_𝑇𝐴/𝜇²𝑦. In the case of 𝑚>1,𝑣 is proportional to 𝑥^(𝑚−1)/3 and is strictly monotone increasing to infinity as 𝑥 tends to infinity, which is not accepted in physics. Therefore, we restrict our investigations for the interval −1<𝑚≤1.

We note that the special case 𝑚=1 do admit explicit solution. In [21, 22] the solution to (2.8), (2.9) is given by𝑓(𝜂)=1−𝑒^−𝜂,(2.13) and easy computation shows thatΘ(𝜂)=Φ(Pr)−Ψ(Pr)𝑒^−𝜂+Ω(Pr)𝑒^−2𝜂,(2.14) with Φ(Pr)=[1/[(Pr−1)/Pr+Pr/(Pr−2)−2]][(Pr−1)/Pr], Ψ(Pr)=2/[(Pr−1)/Pr+Pr/(Pr−2)−2], Ω(Pr)=[1/[(Pr−1)/Pr+Pr/(Pr−2)−2]][Pr/(Pr−2)] is the solution to (2.10), (2.11).

Due to the inherent complexity of such flows, to give exact analytical solutions of Marangoni flows are almost impossible. Exact analytical solutions were given by Magyari and Chamkha for thermosolutal Marangoni convection when the wall temperature and concentration variations are quadratic functions of the location [6].

Our goal is to present approximate exponential series solution to the nonlinear boundary value problem (2.8), (2.9), moreover to (2.10), (2.11). Several values of the power law exponent and Prandtl number are considered. The influences of the effects of these parameters are illustrated.

3. Exponential Series Solution

First, our aim is to determine the approximate local solution of 𝑓(𝜂) to (2.8), (2.9). We replace the condition at infinity by one at 𝜂=0. Therefore, (2.8), (2.9) is converted into an initial value problem of (2.8) with initial conditions𝑓(0)=0,𝑓^′(0)=𝜁,𝑓^(0)=−1.(3.1)

In view of the third of the boundary conditions (2.9), let us take the solution of the initial value problem (2.8), (3.1) in the form𝑓(𝜂)=𝛼𝐴₀+^∞_𝑖=1𝐴_𝑖𝑑^𝑖𝑒^{−𝛼𝜂𝑖},(3.2) where 𝛼>0, 𝐴₀=3/(𝑚+2), 𝐴_𝑖(𝑖=1,2,…) are coefficients and 𝛼>0 and 𝑑 are constants. Conditions in (2.9) yield the following equations:𝛼𝐴₀+^∞_𝑖=1𝐴_𝑖𝑑^𝑖𝛼=0,_3∞_𝑖=1𝑖²𝐴_𝑖𝑑^𝑖=−1.(3.3) It may be remarked that the classic Briot-Bouquet theorem [23] guarantees the existence of formal solutions (3.2) to (2.8), (3.1), the value of 𝐴₀, and also the convergence of formal solutions.

Let us introduce the new variable 𝑍 such as𝑍=𝑑𝑒^−𝛼𝜂.(3.4)

It is evident that the third boundary condition in (2.9) is automatically satisfied. From differential equation (2.8) with (3.2) we get−^∞_𝑖=1𝑖³𝐴_𝑖𝑍^𝑖+𝑚+23𝐴⁰+^∞_𝑖=1𝐴_𝑖𝑍^𝑖^∞_𝑖=1𝑖²𝐴_𝑖𝑍^𝑖−2𝑚+13^∞_𝑖=1𝑖𝐴_𝑖𝑍^𝑖²=0.(3.5)

Equating the coefficients of like powers of 𝑍 one can obtain the expressions for coefficients 𝐴₂,𝐴₃,… with 𝑚 and 𝐴₁:𝐴₂=−112𝐴²¹𝐴(𝑚−1),₃=1216𝐴³¹𝐴(𝑚−1)(𝑚−2),₄=−115552𝐴⁴¹𝐴(𝑚−1)4𝑚²−15𝑚+17₅,=14665600𝐴⁵¹𝐴(𝑚−1)62𝑚³−371𝑚²+757𝑚−610₆,=−146656000𝐴⁶¹𝐴(𝑚−1)32𝑚⁴−257𝑚³+810𝑚²−1171𝑚+730₇,=1740710656000𝐴⁷¹×25742𝑚₅(𝑚−1),⋮−263609𝑚⁴+1108202𝑚³−2419211𝑚²+2737856𝑚−1383380,(3.6) From system (3.3) with the choice of 𝐴₁=1 the parameter values of 𝑑 and 𝛼 can be numerically determined. By these parameters the complete series solution (3.2) is reached. The calcutaled values of 𝑑, 𝛼 and 𝑓^(0) are shown in Table 1, and the variation of 𝑓^(0) with 𝑚 on Figure 1.

Table 1 

Figure 1 

Variation of 𝜁 with 𝑚.

The series forms for 𝑓(𝜂) and 𝑓^(𝜂) are given below for some special values of the exponent 𝑚(𝑚=−0.5;𝑚=0;𝑚=1):𝑓(𝜂𝑚)==2−.01.257∶758380−2.806255505𝑒_{−1.06387919𝜂}+0.9252777520𝑒^{−1.06387919𝜂}²−0.3389803468𝑒^{−1.06387919𝜂}³+0.1266709685𝑒^{−1.06387919𝜂}⁴−0.04756400678𝑒^{−1.06387919𝜂}⁵+0.01788183745𝑒^{−1.06387919𝜂}⁶−0.006724299058𝑒^{−1.06387919𝜂}⁷+0.002528494356𝑒^{−1.06387919𝜂}⁸𝑓_−0.0009506669013𝑒^{−1.06387919𝜂}⁹+0.0003573924970𝑒^{−1.06387919𝜂}¹⁰,(𝜂)=2.985516834𝑒^{−1.06387919𝜂}−1.968767491𝑒^{−1.06387919𝜂}²+1.081902410𝑒^{−1.06387919𝜂}−0.5390504296𝑒^{−1.06387919𝜂}⁴+0.2530117851𝑒^{−1.06387919𝜂}⁵−0.1141446885𝑒^{−1.06387919𝜂}⁶+0.05007689285𝑒^{−1.06387919𝜂}⁷−0.02152010022𝑒^{−1.06387919𝜂}⁸𝑓(𝜂𝑚)==+10.04∶.80109901902245552−519.868𝑒3209870^{−1.06387919𝜂}𝑒⁹−0.0038_{−0.9879394966𝜂}02224402𝑒^{−1.06387919𝜂}¹⁰,+0.2389818639𝑒^{−0.9879394966𝜂}²−0.04524080827𝑒^{−0.9879394966𝜂}³+0.009099650359𝑒^{−0.9879394966𝜂}⁴−0.001854352461𝑒^{−0.9879394966𝜂}⁵+0.0003780882937𝑒^{−0.9879394966𝜂}⁶−0.00007689156390𝑒^{−0.9879394966𝜂}⁷+0.00001559224546𝑒^{−0.9879394966𝜂}⁸𝑓_−0.000003154096351𝑒^{−0.9879394966𝜂}⁹+0.0000006368276490𝑒^{−0.9879394966𝜂}¹⁰,(𝜂)=1.662909512𝑒^{−0.9879394966𝜂}−0.4721992446𝑒^{−0.9879394966𝜂}²+0.1340855441𝑒^{−0.9879394966𝜂}³−0.03595961598𝑒^{−0.9879394966𝜂}⁴+0.009159940185𝑒^{−0.9879394966𝜂}⁵−0.002241170151𝑒^{−0.9879394966𝜂}⁶+0.0005317494905𝑒^{−0.9879394966𝜂}⁷−0.0001232335611𝑒^{−0.9879394966𝜂}⁸𝑓(𝜂𝑚)==+110−∶.0𝑒00_−𝜂02804450725𝑒^{−0.9879394966𝜂}⁹−0.000006291471870𝑒^{−0.9879394966𝜂}¹⁰,𝑓_,(𝜂)=𝑒^−𝜂.(3.7) It can be seen that for the case 𝑚=1 the obtained solution coincides with the exact solution (2.13). The effect of the exponent 𝑚 on the velocity profiles 𝑓^(𝜂) is illustrated in Figure 2. The values of 𝑓^(0)=𝜁 decrease as 𝑚 is changing from negative values to positive ones.

Figure 2 

Variation of 𝑓 with 𝜂.

Applying the series solution for 𝑓 the second-order linear differential equation (2.10) for Θ can be solved similarly, which presents the temperature distribution. Here we define Θ(𝜂) as the seriesΘ(𝜂)=𝐵₀+^∞_𝑖=1𝐵_𝑖𝑑^𝑖𝑒^{−𝛼𝜂𝑖},(3.8) with coefficients 𝐵_𝑖(𝑖=0,1,2,…) and hence the individual coefficients will be determined from differential equation (2.10) with (3.2) as follows𝐵₁=𝐴₁𝐵₀Pr𝐵Pr−1(𝑚+1),₂=112𝐴²¹𝐵₀Pr(Pr−1)(Pr−2)3𝑚²𝐵Pr+𝑚²+6𝑚Pr+3Pr−1₃,=−1216𝐴³¹𝐵₀𝑃𝑟𝐹(Pr,𝑚)=𝑚₃(Pr−1)(Pr−2)(Pr−3)𝐹(Pr,𝑚),⋮−𝑚3Pr²−19Pr−2+𝑚²−14Pr²−20Pr+4,(3.9)

Remark that these coefficients as expressions of 𝐵₀ can be calculated only for noninteger values of the low Prandtl numbers. In (2.11) the second boundary condition is automatically satisfied, and from the first condition coefficient 𝐵₀ is to be determined, that is, from the equation𝐵₀+𝐵₁𝑑+𝐵₂𝑑²+𝐵₃𝑑³+⋯=1,(3.10) together with (3.9) (see Table 2).

Table 2 

For Θ(𝜂) with Prandtl number Pr =298 and three values of 𝑚(𝑚=−0.5;𝑚=0;𝑚=1) the first ten terms are given belowΘ(𝜂𝑚)==−−05.251.∶5866457+690.2257842𝑒_{−1.06387919𝜂}−228.3499003𝑒^{−1.06387919𝜂}²+83.59931856𝑒^{−1.06387919𝜂}³−31.32427417𝑒^{−1.06387919𝜂}⁴+11.72818950𝑒^{−1.06387919𝜂}⁵−4.437848713𝑒^{−1.06387919𝜂}⁶+1.649735001𝑒^{−1.06387919𝜂}⁷−.6357069408𝑒^{−1.06387919𝜂}⁸Θ(𝜂𝑚)==+0−2.∶25236.88147023280461+𝑒433.8991^{−1.06387919𝜂}70⁹4−𝑒0.095491457_{−0.9879394966𝜂}98𝑒^{−1.06387919𝜂}¹⁰,−185.8554249𝑒^{−0.9879394966𝜂}²+23.32308704𝑒^{−0.9879394966𝜂}³−11.64354979𝑒^{−0.9879394966𝜂}⁴−.6456997006𝑒^{−0.9879394966𝜂}⁵−1.774772066𝑒^{−0.9879394966𝜂}⁶−.8708402162𝑒^{−0.9879394966𝜂}⁷−.6992264953𝑒^{−0.9879394966𝜂}⁸Θ(𝜂𝑚)==−173.∶521609496856𝑒^{−0.9879394966𝜂}⁹−.4045737228𝑒^{−0.9879394966𝜂}¹⁰,25−4471504(𝑒^−𝜂)+2220175(𝑒^−𝜂)².(3.11)

It may be noted that the Prandtl number Pr=298 corresponds to the power transformer oil. We point out that for the case 𝑚=1 the solution Θ(𝜂) coincides with the exact solution (2.14).

The effects of the power law exponent 𝑚 and the Prandtl number are exhibited in Figures 3–14, where Pr =0.27 corresponds to the mercury and Pr=0.7 corresponds to the air. Figures 3–5 illustrate the influence of the Prandtl number on the temperature Θ for 𝑚=1. It can be observed in Figure 3 that for low Prandtl numbers 0.27≤Pr≤1.00001 the maximum value of Θ decreases as Pr increases and for a high Prandtl numbers 2.5≤Pr≤7.00001 and 70≤Pr≤298 the maximum value of Θ increases as Pr increases. In all three cases the boundary layer thickness increases as Pr increases. Figures 6–8 depict the effect power exponent 𝑚 for fixed values of Pr. It can be observed in Figures 6 and 7 that the boundary layer thickness increases as 𝑚 increases and the maximum value of Θ decreases as 𝑚 increases for Pr=0.27,2.2, while for high Prandtl number (Pr=298) the reverse effect of 𝑚 on the maximum of Θ can be seen. Figures 9–11 illustrate the effect of Pr on Θ^′ for 𝑚=1 and Figures 12–14 represent the effect of 𝑚 for Θ^ for Pr=0.27,2.2,298.

Figure 3 

Variation of Θ with Pr (0.27≤Pr≤1.00001) for 𝑚=1.

Figure 4 

Variation of Θ with Pr (2.5≤Pr≤7.00001) for 𝑚=1.

Figure 5 

Variation of Θ with Pr (70≤Pr≤298) for 𝑚=1.

Figure 6 

Variation of Θ for Pr=0.27.

Figure 7 

Variation of Θ for Pr=2.2.

Figure 8 

Variation of Θ for Pr=298.

Figure 9 

The effect of the Prandtl number on Θ′ for 𝑚=1 (0.27≤Pr≤1.00001).

Figure 10 

The effect of the Prandtl number on Θ′ for 𝑚=1 (2.5≤Pr≤7.00001).

Figure 11 

The effect of the Prandtl number on Θ′ for 𝑚=1 (70≤Pr≤298).

Figure 12 

The effect of 𝑚 on Θ′ for Pr=0.27 (𝑚=−0.5;0;1).

Figure 13 

The effect of 𝑚 on Θ′ for Pr=2.2 (𝑚=−0.5;0;1).

Figure 14 

The effect of 𝑚 on Θ′ for Pr=298 (𝑚=−0.5;0;1).

4. Conclusion

In this paper the incompressible flow and heat transfer over a flat impermeable plate has been investigated. The resulting governing equations have been transformed into a system of nonlinear ordinary differential equations by applying suitable similarity transformation. For these equations approximate exponential series solutions are determined and the effects of the power exponent and the Prandtl number are illustrated in Figures 1–14, and it is observed that the values of 𝑓^′ decrease as power exponent 𝑚 increases (see Figures 1 and 2), moreover the boundary layer thickness increases as 𝑚 or Pr increases. The temperature profiles are exhibited in Figures 3–8, it is observed that for low Prandtl number the temperature Θ decreases as Pr increases and for high Prandtl numbers the influence of Pr is opposite.

Acknowledgment

This research was carried out as part of the TAMOP-4.2.1.B-10/2/KONV-2010-0001 Project with support by the European Union, cofinanced by the European Social Fund.