2.1. Governing equations and boundary conditions

Figure 1 schematically illustrates the problem under present investigation. We study the Marangoni instability in a thin layer of an incompressible viscoelastic polymer solution, initially resting on a flat rigid substrate (of lower thermal conductivity compared with the liquid) in the gravitational field g. This laterally infinite, two-dimensional film is separated from the ambient gas phase by its deformable free surface located at $z = h(x,t)$. The polymeric solution is defined by its relaxation time $\rlap{-}{\lambda }$, viscosity ${\mu _o}$, density $\rho$, thermal conductivity $\kappa$, thermal diffusivity $\alpha$, mass diffusivity D and surface tension $\sigma$.

Figure 1.

Schematic of the physical system under investigation. Marangoni instability is induced in a thin viscoelastic polymer film confined between its deformable free surface (located at $z = h(x,t)$) and a flat substrate (at the $z = 0$ plane) when subjected to heating from above. The polymeric solution is a binary mixture of a Newtonian solvent with a polymeric solute. The incorporation of the Soret effect signifies the combined thermosolutal instability in the system.

The entire liquid film is subjected to a uniform transverse temperature gradient, specified to be $- \vartheta$ at the $z = 0$ plane. Thus, a negative (positive) $\vartheta$ indicates the case of heating the film from the air–liquid interface (substrate). Here, we are interested in investigating the instability phenomenon only for the former case (i.e. heating the fluid layer at the air–liquid interface). This applied temperature gradient induces a concentration gradient in the film via the Soret effect. The heat and mass fluxes in the bulk of the liquid layer are thus given by (de Groot & Mazur Reference De Groot and Mazur2011)

(2.1)\begin{gather}{\boldsymbol{J}_H} = - \kappa {\kern 1pt} {\boldsymbol{\nabla} }T,\end{gather}

(2.2)\begin{gather}{\boldsymbol{J}_M} = - \rho D({\boldsymbol{\nabla} }c + {{\mathcal S}}{\kern 1pt} {\boldsymbol{\nabla} }T),\end{gather}

respectively, where T is the temperature, c is the solute concentration and ${{\mathcal S}}$ is the Soret diffusion coefficient of the mixture. For a polymeric solution, ${{\mathcal S}}$ can be either positive or negative depending upon the solvent quality and the mole fractions of the components, as well as being based on the temperature of the mixture, as mentioned in § 1. However, the Dufour effect, through which the concentration gradient couples back to the dynamics of the temperature field, is neglected in this analysis owing to its exceedingly weak impact on liquids. Equations (2.1) and (2.2) indicate that, in the conductive state with H as the unperturbed film thickness, the applied heat flux generates a temperature difference $\Delta T = |\vartheta |H$ across the film, which, in turn, produces a concentration difference $c = - {{\mathcal S}}\Delta T$ via the Soret effect.

Now, above a certain critical temperature gradient, the thermocapillary and solutocapillary forces on the free liquid surface induce Marangoni convection in the liquid film. Here, we assume the surface tension to vary monotonically with temperature and concentration of the mixture, dictated by the relationship

(2.3)\begin{equation}\sigma = {\sigma _o} - {\sigma _T}(T - {T_r}) + {\sigma _c}(c - {c_r}),\end{equation}

where ${\sigma _o}$ is the surface tension at the reference temperature ${T_r}$ and concentration ${c_r}$. For most polymer blends ${\sigma _T}\;( = - \partial \sigma /\partial T{|_{T = {T_r}}}) \gt 0$, while ${\sigma _c}\;( = \partial \sigma /\partial c{|_{c = {c_r}}})$ can be either positive or negative. The effects of buoyancy are neglected in this study considering the small thickness of the film ($H\lesssim O(1)\;\textrm{cm}$; Pearson Reference Pearson1958). Furthermore, except for $\sigma$, all other thermophysical properties are assumed to remain invariant throughout the analysis. Therefore, the evolutions of the film velocity $\boldsymbol{v} \equiv \{ u(x,z,t),w(x,z,t)\}$, pressure $p(x,z,t)$, temperature $T(x,z,t)$ and solute concentration $c(x,z,t)$ with time t over the horizontal range $x \in ( - \infty ,\infty )$ and the vertical range $z \in [0,h]$ are governed by

(2.4)\begin{gather}{\boldsymbol{\nabla\, \cdot\, }}{\kern 1pt} \boldsymbol{v} = \boldsymbol{0},\end{gather}

(2.5)\begin{gather}\rho \left( {\frac{{\partial \boldsymbol{v}}}{{\partial t}} + \boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }}\boldsymbol{v}} \right) = - {\boldsymbol{\nabla} }p + {\boldsymbol{\nabla} \,\boldsymbol{\cdot}\, }\boldsymbol{\tau } - \rho g \boldsymbol{k},\end{gather}

(2.6)\begin{gather}\frac{{\partial T}}{{\partial t}} + \boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }}T = \alpha {\nabla ^2}T,\end{gather}

(2.7)\begin{gather}\frac{{\partial c}}{{\partial t}} + \boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }}c = D{\nabla ^2}c + {{\mathcal S}}D{\nabla ^2}T,\end{gather}

respectively, where $\boldsymbol{\tau } = \left[ {\begin{aligned} {{\boldsymbol{\tau }_{xx}}}&\quad {{\boldsymbol{\tau }_{xz}}}\\ {{\boldsymbol{\tau }_{zx}}}&\quad {{\boldsymbol{\tau }_{zz}}} \end{aligned}} \right]$ is the deviatoric stress tensor, k is the unit vector in the z-direction and ${\boldsymbol{\nabla} } \equiv \{ \partial /\partial x,\partial /\partial z\}$. Note that the dynamics of the gas and liquid phases are decoupled here considering the large ratios between their densities, viscosities and thermal diffusivities.

The boundary conditions that accompany the set of governing equations (2.4)–(2.7) are as follows. At the rigid substrate, we impose the no-slip, no-penetration condition for velocity, a specified heat flux and the mass impermeability conditions, represented by

(2.8a–c)\begin{equation}\boldsymbol{v} = \boldsymbol{0},\quad \frac{{\partial T}}{{\partial z}} = - \vartheta ,\quad \frac{{\partial c}}{{\partial z}} = {{\mathcal S}}\vartheta \quad \textrm{at}\;z = 0,\end{equation}

respectively.

At the deformable free surface, i.e. at $z = h(x,t)$, the boundary conditions comprise the kinematic condition

(2.9a)\begin{equation}w = \frac{{\partial h}}{{\partial t}} + u\frac{{\partial h}}{{\partial x}},\end{equation}

which states that the velocity of the free surface is equal to the velocity of the liquid, thus giving its location. The balance of the tangential and normal stress components at the free surface reads

(2.9b)\begin{gather}\frac{1}{{\sqrt {1 + {{(\partial h/\partial x)}^2}} }}\left\{ {{\tau_{xz}}\left[ {1 - {{\left( {\frac{{\partial h}}{{\partial x}}} \right)}^2}} \right] + {\tau_{zz}}\frac{{\partial h}}{{\partial x}} - {\tau_{xx}}\frac{{\partial h}}{{\partial x}}} \right\} = \frac{{\partial \sigma }}{{\partial x}} + \frac{{\partial \sigma }}{{\partial z}}\frac{{\partial h}}{{\partial x}},\end{gather}

(2.9c)\begin{gather}- p + \frac{1}{{1 + {{(\partial h/\partial x)}^2}}}\left[ {{\tau_{zz}} + {\tau_{xx}}{{\left( {\frac{{\partial h}}{{\partial x}}} \right)}^2} - 2{\tau_{xz}}\frac{{\partial h}}{{\partial x}}} \right] = \sigma {\kern 1pt} {{\mathcal H}},\end{gather}

where ${{\mathcal H}} = ({\partial ^2}h/\partial {x^2}){[1 + {(\partial h/\partial x)^2}]^{ - 3/2}}$ is the mean curvature.

The thermal boundary condition at the free surface includes the balancing of heat flux across the interface. This heat exchange process with the ambient gas phase is approximated in this analysis by the heat transfer coefficient q between the liquid and the gas phase as follows:

(2.9d)\begin{equation}- \kappa \left( {\frac{{\partial h}}{{\partial x}}\frac{{\partial T}}{{\partial x}} - \frac{{\partial T}}{{\partial z}}} \right) + q(T - {T_\infty })\sqrt {1 + {{(\partial h/\partial x)}^2}} = 0,\end{equation}

where ${T_\infty }$ is the uniform gas temperature.

Finally, for this non-volatile binary mixture, the mass flux vanishes at the free surface. Mathematically, this is expressed by

(2.9e)\begin{equation}\kappa \left( { - \frac{{\partial h}}{{\partial x}}\frac{{\partial c}}{{\partial x}} + \frac{{\partial c}}{{\partial z}}} \right) - {{\mathcal S}}q(T - {T_\infty })\sqrt {1 + {{(\partial h/\partial x)}^2}} = 0.\end{equation}

2.2. Constitutive equation for the fluid

Viscoelastic fluids exhibit complex rheology owing to both the viscous and elastic properties. To depict the rheology of these fluids, a wide variety of constitutive models have been developed over the years. The rheology of the fluid is approximated in this work by the Maxwell model (Bird et al. Reference Bird, Armstrong and Hassager1987)

(2.10)\begin{equation}\boldsymbol{\tau } + \rlap{-}{\lambda } \frac{{\textrm{D}\boldsymbol{\tau }}}{{\textrm{D}t}} = {\mu _o}[({\boldsymbol{\nabla} }\boldsymbol{v}) + {({\boldsymbol{\nabla} }\boldsymbol{v})^\textrm{T}}],\end{equation}

which characterizes the fluid by a single relaxation time $\rlap{-}{\lambda } \;( = {\mu _o}/G$, where G is the elastic modulus). Out of the spectrum of relaxation time exhibited by the liquid, $\rlap{-}{\lambda }$ in (2.10) is interpreted as the longest relaxation time. Moreover, in (2.10),

(2.11)\begin{equation}\frac{{\textrm{D}\boldsymbol{\tau }}}{{\textrm{D}t}} = \frac{{\partial \boldsymbol{\tau }}}{{\partial t}} + (\boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }})\boldsymbol{\tau } - {({\boldsymbol{\nabla} }\boldsymbol{v})^\textrm{T}}\,\boldsymbol{\cdot }\,\boldsymbol{\tau } - \boldsymbol{\tau }\,\boldsymbol{\cdot} \,({\boldsymbol{\nabla} }\boldsymbol{v})\end{equation}

is the upper convected derivative. However, in this investigation, since a linear stability analysis will be carried out for small perturbations around an initially quiescent state, the nonlinear terms in (2.11) will not make any contributions towards the final results. Hence, $\textrm{D}\boldsymbol{\tau }/\textrm{D}t \equiv \partial \boldsymbol{\tau }/\partial t$ for the present analysis. From (2.10) it is clear that, in the limit $\rlap{-}{\lambda } \to 0$, the Maxwell model depicts the Newtonian fluid behaviour.

2.3. Non-dimensionalization

Let us now non-dimensionalize the boundary value problem (BVP) formulated by (2.4)–(2.9). Considering the unperturbed film thickness H as the characteristic length scale, the thermal diffusion time ${H^2}/\alpha$ as the characteristic time scale and $|\vartheta |H$ as the temperature scale, we define the following set of dimensionless variables:

(2.12) \begin{equation}\left. \begin{array}{l} (\bar{x},\bar{z}) = \dfrac{{(x,z)}}{H},\quad \bar{h} = \dfrac{h}{H},\quad \bar{t} = \dfrac{t}{{{H^2}/\alpha }},\quad (\bar{u},\bar{w}) = \dfrac{{u,w}}{{(\alpha /H)}},\quad \bar{\tau } = \dfrac{\tau }{{\mu_{{o}} \alpha /{H^2}}},\\ \bar{p} = \dfrac{p}{{\mu_{{o}} \alpha /{H^2}}},\quad \bar{T} = \dfrac{{T - {T_\infty }}}{{|\vartheta |H}},\quad \bar{c} = \dfrac{c}{{{\sigma_T}|\vartheta |H/{\sigma_c}}}. \end{array} \right\}\end{equation}

It may be noted that the characteristic scale adopted in (2.12) coincides with the previous works by Pearson (Reference Pearson1958), Shklyaev, Nepomnyashchy & Oron (Reference Shklyaev, Nepomnyashchy and Oron2009), helping to compare this work with these previously reported studies. Dropping the overbar from the non-dimensional variables for convenience in presentation, we finally arrive at the following set of dimensionless governing equations:

(2.13)\begin{gather}{\boldsymbol{\nabla}\, \boldsymbol{\cdot}\, }\boldsymbol{v} = \boldsymbol{0},\end{gather}

(2.14)\begin{gather}P{r^{ - 1}}\left( {\frac{{\partial {\kern 1pt} \boldsymbol{v}}}{{\partial {\kern 1pt} t}} + \boldsymbol{v}\,\boldsymbol{\cdot\, \nabla }\boldsymbol{v}} \right) = - {\boldsymbol{\nabla} }p + {\boldsymbol{\nabla}\,\boldsymbol{\cdot}\, }\boldsymbol{\tau } - Ga\,\boldsymbol{k},\end{gather}

(2.15)\begin{gather}\frac{{\partial T}}{{\partial t}} + \boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }}T = {\nabla ^2}T,\end{gather}

(2.16)\begin{gather}\frac{{\partial c}}{{\partial {\kern 1pt} t}} + \boldsymbol{v}{\boldsymbol{\,\cdot\, \nabla }}c = Le({\nabla ^2}c + \chi {\nabla ^2}T).\end{gather}

The boundary conditions (2.8) and (2.9) now take the form

(2.17a–c)\begin{equation}\boldsymbol{v} = \boldsymbol{0},\quad \frac{{\partial T}}{{\partial z}} = - {{\mathcal Q}},\quad \frac{{\partial c}}{{\partial z}} = \chi {{\mathcal Q}}\quad \textrm{at}\;z = 0,\end{equation}

and

(2.18a)\begin{gather}w = \frac{{\partial h}}{{\partial t}} + u\frac{{\partial h}}{{\partial x}},\end{gather}

(2.18b)\begin{gather}\left( {\frac{{\partial T}}{{\partial z}} - \frac{{\partial h}}{{\partial x}}\frac{{\partial T}}{{\partial x}}} \right) + Bi\,T\sqrt {1 + {{(\partial h/\partial x)}^2}} = 0,\end{gather}

(2.18c)\begin{gather}\left( {\frac{{\partial c}}{{\partial z}} - \frac{{\partial h}}{{\partial x}}\frac{{\partial c}}{{\partial x}}} \right) - \chi Bi\,T\sqrt {1 + {{(\partial h/\partial x)}^2}} = 0,\end{gather}

(2.18d)\begin{gather}- p + \frac{1}{{1 + {{(\partial h/\partial x)}^2}}}\left[ {{\tau_{zz}} + {\tau_{xx}}{{\left( {\frac{{\partial h}}{{\partial x}}} \right)}^2} - 2{\kern 1pt} {\tau_{xz}}\frac{{\partial h}}{{\partial x}}} \right] = \Sigma \frac{{{\partial ^2}h/\partial {x^2}}}{{{{[1 + {{(\partial h/\partial x)}^2}]}^{3/2}}}},\end{gather}

(2.18e)\begin{gather}\begin{array}{ccccc} & \dfrac{1}{{\sqrt {1 + {{(\partial h/\partial x)}^2}} }}\left\{ {{\tau_{xz}}\left[ {1 - {{\left( {\dfrac{{\partial h}}{{\partial x}}} \right)}^2}} \right] + {\tau_{zz}}\dfrac{{\partial h}}{{\partial x}} - {\tau_{xx}}\dfrac{{\partial h}}{{\partial x}}} \right\}\\ & \quad = Ma\left[ {\dfrac{{\partial {\kern 1pt} c}}{{\partial {\kern 1pt} x}} - \dfrac{{\partial T}}{{\partial x}} + {\kern 1pt} \left( {\dfrac{{\partial {\kern 1pt} c}}{{\partial z}} - \dfrac{{\partial T}}{{\partial z}}} \right)\dfrac{{\partial h}}{{\partial {\kern 1pt} x}}} \right]\quad \textrm{at}\;z = h(x,t), \end{array}\end{gather}

while the Maxwell constitutive model (2.10) reads

(2.19)\begin{equation}\boldsymbol{\tau } + De\frac{{\partial \boldsymbol{\tau }}}{{\partial t}} = ({\boldsymbol{\nabla} }\boldsymbol{v}) + {({\boldsymbol{\nabla} }\boldsymbol{v})^\textrm{T}}.\end{equation}

Note that, the non-dimensional parameter ${{\mathcal Q}} = \vartheta /|\vartheta |$ introduced in (2.17) indicates the direction of applied temperature gradient. The value ${{\mathcal Q}} = \textrm{1}$ represents the case of heating the fluid layer from below, while ${{\mathcal Q}} = - \textrm{1}$ stands for heating from above. Besides ${{\mathcal Q}}$, the BVP (2.13)–(2.19) is further governed by the following set of dimensionless parameters: the Marangoni number Ma, the Soret number $\chi$, the Deborah number De, the (inverse) Lewis number Le, the Biot number Bi, the Prandtl number Pr, the Galileo number Ga, and the (inverse) capillary number $\Sigma$:

(2.20)\begin{equation}\left. \begin{array}{l} Ma = \dfrac{{{\sigma_T}|\vartheta |{H^2}}}{{{\mu_o}\alpha }}\textrm{,}\quad \chi = \dfrac{{{{\mathcal S}}{\sigma_c}}}{{{\sigma_T}}},\quad De = \dfrac{{\rlap{-}{\lambda } \alpha }}{{{H^2}}},\quad Le = \dfrac{D}{\alpha },\\ Bi = \dfrac{{qH}}{\kappa },\quad Pr = \dfrac{{{\mu_o}}}{{\rho \alpha }},\quad Ga = \dfrac{{\rho g{H^3}}}{{{\mu_o}\alpha }},\quad \Sigma = \dfrac{{\sigma H}}{{{\mu_o}\alpha }}. \end{array} \right\}\end{equation}

The Marangoni number governs the present instability phenomenon. It gives the critical temperature difference across the film (|ϑ|H) for which convection sets in, overcoming the stabilizing actions of viscous and thermal diffusion. The Soret number takes into account the relative contributions of the thermocapillary and solutocapillary forces towards the free surface force. Depending upon the Soret coefficient ${{\mathcal S}}$, $\chi$ can assume both positive and negative values for a polymeric mixture. The Deborah number quantifies the elastic behaviour of the fluid through the magnitude of $\rlap{-}{\lambda }$. In this analysis, $De = 0$ indicates a Newtonian binary mixture ($\rlap{-}{\lambda } = 0$), whereas increasing values of De signify enhanced elasticity of the mixture. The (inverse) Lewis number compares the characteristic mass diffusion time scale $({H^2}/D)$ with the thermal diffusion time scale $({H^2}/\alpha )$. The Biot number characterizes the heat transfer rate across the free surface. The Prandtl number is a material property of the fluid, representing the ratio between thermal diffusion time scale $({H^2}/\alpha )$ and viscous diffusion time scale $(\rho {H^2}/{\mu _o})$. The Galileo number and the (inverse) capillary number take into account the deformability of the free surface through the magnitude of g and $\sigma$, respectively.

Formulating the problem, we now conclude this section by discussing the physically permissible range of the above-mentioned non-dimensional parameters. For a viscoelastic binary mixture, $Pr \gg 1$ and Le typically ranges between $O({10^{ - 5}})\lesssim Le\lesssim O({10^{ - 1}})$. Furthermore, we analyse here both the separate cases of $\chi \gt 0$ and $\chi \lt 0$, and vary Bi within the range $0 \lt Bi \lt 1$. To study the stability characteristics of both the weakly and highly viscoelastic fluids, we consider a broad spectrum for De: $0\lesssim De\lesssim O(10)$. Note that, for a 0.1 mm thick polymeric film with $\alpha \approx O({10^{ - 7}})\;{\textrm{m}^2}\;{\textrm{s}^{ - 1}}$, this range of De encompasses the fluids having $\rlap{-}{\lambda } \approx O(0\textrm{--}10)\;\textrm{s}$. Typical values of $\rlap{-}{\lambda }$ for different polymeric solutions are given in table 1.

Table 1.

Typical values of relaxation time for different viscoelastic fluids (from Adam & Delsanti Reference Adam and Delsanti1983; Joseph Reference Joseph1990; Ebagninin, Benchabane & Bekkour Reference Ebagninin, Benchabane and Bekkour2009).

Moreover, to understand the role of surface deformability on the onset of instability in the system, we consider here both the cases of a deformable and a non-deformable free surface. Since increasing gravitational and surface tension forces reduce the deformability of a free surface, we consider the conditions $(Ga,\Sigma ) = (0.1,{10^3})$ to depict a deformable free surface. For a 0.1 mm thick layer of the polymeric solution with $\rho \approx O({10^3})\;\textrm{kg}\;{\textrm{m}^{ - 3}}$ and ${\mu _o} \approx O({10^{ - 2}})\;\textrm{Pa}\;\textrm{s}$, the above-mentioned values of Ga and $\Sigma$ refer to $g \approx 0.1\;\textrm{m}\;{\textrm{s}^{ - 2}}$ and $\sigma \approx 0.01\;\textrm{N}\;{\textrm{m}^{ - 1}}$. On the other hand, the free surface is treated as non-deformable in the limit $(Ga,\Sigma ) \to \infty$, which typically refers to a liquid layer with significantly high surface tension in the terrestrial environment.

4.1. Validation of the numerical scheme

Before proceeding further, let us first verify the accuracy of the employed numerical scheme. Note that, although a few studies on the Marangoni instability in a pure viscoelastic liquid film have been reported in the literature (Getachew & Rosenblat Reference Getachew and Rosenblat1985; Dauby et al. Reference Dauby, Parmentier, Lebon and Grmela1993; Parmentier et al. Reference Parmentier, Lebon and Regnier2000), an investigation considering the binary aspect of the fluid has not been carried out before. Furthermore, in contrast to these works, a slightly different system configuration is adopted in the present work (instead of temperature, its gradient at the solid substrate is specified in this problem). Therefore, we first test the validity of our numerical scheme for a system with boundary conditions identical to the present study. In figure 2(a) the results of the present computations are compared against the well-known results of Pearson (Reference Pearson1958) (for the monotonic mode in purely thermocapillary-driven instability in a Newtonian liquid film, i.e. for $(\varOmega ,De,\chi ,{{\mathcal Q}}) = (0,0,0,1)$) and Shklyaev et al. (Reference Shklyaev, Nepomnyashchy and Oron2009) (for the oscillatory mode in thermosolutal Marangoni instability in a Newtonian binary liquid film, i.e. for $(De,{{\mathcal Q}}) = (0,1)$). It can be observed that an excellent quantitative agreement exists between the results of the present numerical implementation and the above-cited papers for the entire range of wavenumbers.

Figure 2.

Validation of the numerical method. (a) Comparison of the present results with Pearson (Reference Pearson1958) and Shklyaev et al. (Reference Shklyaev, Nepomnyashchy and Oron2009) via the neutral stability curve for the monotonic and oscillatory instability mode. To replicate the case of a Newtonian liquid film with an insulated, non-deformable free surface subjected to heating from below, we consider $(Bi,De) = 0$, $(Ga,\Sigma ) \to \infty$ and ${{\mathcal Q}} = 1$. For the binary mixture: $Pr = 2,$ $\chi = - 0.2$ and $Le = {10^{ - 3}}$. (b) Comparison of the present results with Getachew & Rosenblat (Reference Getachew and Rosenblat1985) for pure thermocapillary-driven instability in a viscoelastic liquid film. To reproduce the results of the above-cited paper, we consider $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \theta } = 0$ in (3.20c), $(Bi,\chi ) = 0$, ${{\mathcal Q}} = 1$ and $(Ga,\Sigma ) \to \infty$. Lines marked by 1, 2, 3 and 4 correspond to $(De,Pr) = (0.12,1)$, $(0.45,\,{\kern 1pt} 1)$, $(0.5,\,{\kern 1pt} 0.1)$ and $(0.1,\,{\kern 1pt} 10)$, respectively.

Another attempt is made in figure 2(b) to verify the accuracy of our numerical scheme, this time for a viscoelastic fluid, against the results of Getachew & Rosenblat (Reference Getachew and Rosenblat1985). In order to do so, we consider $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \theta } = 0$ in (3.20c) and $(\chi ,{{\mathcal Q}}) = (0,1)$. In figure 2(b), we compare our results with the numerical computations of Getachew & Rosenblat (Reference Getachew and Rosenblat1985) for a wide range of $De\;( = (0.12,\,{\kern 1pt} 0.45)$, top panel) and $Pr\;( = (0.1,\,{\kern 1pt} 10)$, bottom panel). As can be seen from figure 2(b), a fairly good agreement exists between the results for each value of the parameter set (De, Pr). These comparisons ensure the accuracy of the present numerical implementation.