# Influence of Temperature Gradient on the Stability of a Viscoelastic Film with Gradually Fading Memory over a Topography with Ridges and Furrows

Authors: Mohammed Rizwan Sadiq Iqbal
Journal: Chinese Journal of Engineering
 (2013)

---

## Abstract

A long-wave evolution equation is derived using an asymptotic analysis, and the linear stability of a viscoelastic film flowing along the direction of parallel grooves over a uniformly heated topography is studied. A numerical approach adopting spectral collocation technique is used to demarcate the stable and unstable flow regimes. The combined influence of thermocapillarity and viscoelasticity on the films stability is analyzed. By accounting the bottom topography comprising longitudinal gutters, the possibilities of regulating the film dynamics under iso- and nonisothermal conditions and/or optimizing design structure of an apparatus for desirable flow outcomes have been focussed.


---

## Body

## 1. Introduction

Falling films are relevant to a broad class of interfacial instability problems over a wide range of length and time scales in various technological setups since they offer small thermal resistance and large contact area at small specific flow rates. This serves as an advantage in many processes involving cooling, condensation, absorption, and evaporation, thereby regulating the transport of mass, momentum, and heat across the liquid-gas and liquid-solid interfaces. For a detailed review of falling film problems, refer to Nepomnyashchy et al. [1], Kalliadasis et al. [2], and the references therein.Although many studies available in literature have focused attention on Newtonian fluids, many fluids used in industrial applications are rheologically complex and non-Newtonian in nature. This has led to the generalization of Navier-Stokes model to satisfy highly nonlinear constitutive laws to arrive at complex partial differential equations, which are one order higher than the Navier-Stokes equations [3–8].Unlike Newtonian fluids, which respond virtually instantaneously to an imposed deformation rate, viscoelastic fluids respond on a macroscopically large time scale, known as the relaxation time. The viscoelastic fluids, a subclass of microstructure flows, display both elastic (for deformation rates larger than the inverse relaxation time) and viscous (for deformation rates smaller than the magnitude of the inverse relaxation time) characteristics. The stress in this liquid is neither directly proportional to the strain nor to the rate of strain but displays a complex relationship [9]. A lot of work on the flow and heat transfer characteristics of non-Newtonian fluids has also been done in order to control the quality of the end product in many manufacturing and processing industries ([10] and the references therein). This area of fluid dynamics can simulate accurately many complex polymeric, biotechnological, and tribological fluids and has attracted special attention with the advent and growth of synthetic material industries. It finds applications in clothing, food, and pharmaceutical industry, and in analyzing the lubrication behavior of bearings, gears, and cams.There are various constitutive models which address the elasticoviscous aspect of a viscoelastic fluid. The Walters-B liquid model [11] can simulate accurately many complex polymeric, biotechnological, and tribological fluids and forms an important basis for the manufacture of plastic engineering equipments, contact lenses, and so forth. The mixture of polymethyl methacrylate and pyridine at 25°C containing 30.5 g of polymer per liter behaves very close to this model. The physical properties of such a fluid constitute its density ρ=0.98×103Kg/m3, limiting viscosity μ=0.79Ns/m2, kinematic viscosity ν=806×10-6Nsm/Kg, surface tension σ=40×10-3N/m, and viscoelastic coefficient Γ0=0.04Ns2/m2. Hydrodynamic stability studies for the mixture of above polymer have been reported by Andersson and Dahl [12], Sadiq and Usha [13], Lin and Chen [14], and Sirwah and Zakaria [15] for thin films.Chang [16] in his report categorized different types of waves while observing a falling film. He identified interesting regimes where the wave behavior changes depending upon the magnitude of Reynolds number, Re. In the low-moderate Reynolds number regime, long gravity-capillary waves dominate. This region is interesting both from a practical point of view (since many industrial equipments operate under this condition) as well as from a theoretical point of consideration (due to regular structures, it is possible to describe the aspects of the flow by using relatively simple models).Tackling the complete set of mass, momentum, and energy equations poses a formidable challenge. Therefore, with lubrication approximation at the forefront, dimension reduction techniques that take advantage of thin film geometry have been successfully developed in the majority of theoretical work that has appeared to date. This assumption neglects relatively unimportant terms while scaling complex equations by considering the depth-length ratio to be small and is an essential analytical tool in modeling thin film low-Reynolds number flows (a subregime of1<Re<300). Under this assumption, the physical variables are expanded in terms of the small aspect ratio to yield a more tractable yet highly nonlinear degenerate fourth-order PDE for the local layer, which is still capable of retaining the dominant physics [17]. For nonisothermal problems, some notable pieces of research include Joo et al. [18], Miladinova et al. [19], Scheid [20], and Sadiq et al. [21].Fluid actuation and transport are core functionalities and challenges in many fluidic systems. The shape and property of the substrate (its curvature and/or flexibility, patterns and/or obstacles) play an important role in the development of flow instabilities and affect the film dynamics and heat transfer. In particular, the interplay between the physical surface and the local fluid flow affects the entire system, where the viscosity, surface tension, thermocapillarity, and the local interface, for instance decides the flow mechanism. A better understanding of these properties and how they can be harnessed and controlled is extremely important. Their impact in relation to optimizing resources and manufacturing cost in the industrial sector have yet to be fully realized, with their use in relation to forming novel functional coatings for applications in many engineering disciplines still at an early stage.There are several investigations on the flow of a thin filmacross a sinusoidal geometry for isothermal [22–30] and non-isothermal problems [31–34]. Another direction of study includes flows along the grooves of corrugated geometries [35–38]. Recently, Heining et al. [39] reported on the flow over different three-dimensional periodic corrugations and found that the interplay between the geometry of the topography and the inertia of the film exhibits a rich variety of pattern formation phenomena.Literature on non-Newtonian fluids, their stability, and dynamics on substrates with topography is not abundant. For a power-law based non-Newtonian model, the effect of substrate topography, inertia, and non-Newtonian rheology was examined for a two-dimensional thin Carreau-Bird model based liquid which has direct relevance in problems connected to the early stages of coating process [40]. The effect of local microtopography with steps, hills, and periodicity on the steady flow was considered by Mogilevskii and Shkadov [41]. The effect of inertia and surface tension on the flow of a power-law film over an undulated surface has been investigated by Heining and Aksel [30]. They identified that the shear-thinning rheology of the film leads to possible resonance. Usha and Uma [42] extended the study carried by Tougou [43] to viscoelastic Walters-B liquid on a wavy incline and found the impact of the viscoelastic parameter on the films stability and long waves generated. Using a convected Maxwell model, Saprykin et al. [44] investigated the influence of inertia and viscoelasticity on a film flowing on a step-down topography. The effect of viscoelastic properties under creeping flow conditions was investigated by Pavlidis et al. [45] considering a Phan-Thien and Tanner viscoelastic model for fluid flowing on a horizontal trench. By incorporating an asymptotic expansion of the upper-convected Maxwell model for elastic fluids, recently, Bouchut and Boyaval [46] proposed a new model for gravity-driven free-surface flows of shallow viscoelastic fluids across a topography. An important feature of their model is the convexity with respect to the physically relevant pseudoconservative variables.Shaqfeh et al. [47] assessed the stability of gravity driven viscoelastic flow of an Oldroyd-B model in a broad spectrum of dimensionless groups categorizing the flow into thick dilute, thin dilute and thick concentrated or melt film regimes. In the present problem, the influence of non-isothermal effects and bottom topography is considered on a viscoelastic fluid’s stability. For this purpose, a Walters-B liquid is considered to model viscoelastic effects. Though this investigation does not cover the entire class of viscoelastic flows (thin to thick films), yet it gives a qualitative understanding and a first insight into the mechanism occurring in thin film regime, where the surface-tension forces become very important. Also, it forms as the first study in modeling such fluids on a corrugated geometry comprising longitudinal grooves. To cater non-isothermal effects, a uniform temperature (larger than the constant ambient air temperature) is prescribed on the bottom of the topography, and a long-wave equation is derived. The stability of the primary flow to infinitesimal disturbances is discussed under non-isothermal effects.


## 2. Problem Description

Figure1 displays the geometry of the physical problem under consideration. It is assumed that the flow rate is high so that the topography is inundated with nonvolatile fluid. The flow is oriented along the grooves in the positive x-direction, and the geometry extends infinitely along the x- and y-directions. The bottom topography, l(x,y)=alcos(2πy/λ) undulates weakly along the y-axis such that f(x,y,t)=l(x,y)+h(x,y,t) measures the total thickness. The parameters al and λ represent the amplitude (the peak deviation of l from its center position) and the wavelength of a typical periodic groove. The air above the liquid is passive such that its temperature, T∞, and pressure, p∞, are constants sufficiently far from the interface. A uniform temperature distribution TW is prescribed at the bottom of the topography assuming the substrate to be a perfect conductor of heat of infinite heat capacity. Boussinesq approximation is used to describe the bulk movement of the liquid. In the framework of this approximation, density (ρ), kinematic viscosity (ν), thermal conductivity (λc) and thermal diffusivity (χ) are assumed constants [1, 20, 48]. The surface tension decreases linearly as a function of temperature and is given by σ = σ∞-γ(T-T∞), where σ∞ is the mean surface tension at gas temperature and γ=-dσ/dT is a positive constant as for most common liquids.Figure 1
Schematic representation of the flow of a viscoelastic film over a topography comprising furrows and ridges.Constitutive equations which result from a retarded-motion expansion represent the correct constitutive equations for flows in which the rate-of-strain tensor and its time derivatives are small [49]. The second order fluid is derived through a retarded motion expansion and shows two different normal stress functions [49] in simple shear flows, being valid for slow and weak flows [50]. Though the shear thinning seen in experiments on polymeric fluids is not described by this model, yet it finds usage in many applications, and it also predicts most non-Newtonian effects qualitatively, if not always quantitatively [51].Dandapat and Samanta [52] assessed the constitutive equations of the second order incompressible fluid [49] framed on the postulate of gradually fading memory [53]. Such fluids show normal stress differences in the shear flow (which is a characteristic property of a viscoelastic fluid), and the constitutive equation is applicable to the flow of dilute polymer solutions (e.g., polyethylene oxide in aqueous solution (POLYOX), polyisobutylene in cetane), which are only slightly viscoelastic and may have a bearing on the study of fluids with drag-reducing properties (Toms phenomenon). Regarding the applicability of their model, if the disturbance time scale is large compared with the characteristic time scale (relaxation time) of the fluid, then the second order fluid model is internally consistent with the stress-relaxing fluid due to Oldroyd [3]. As pointed out by Porteous and Denn [54], this would happen if the viscoelastic parameter is small.The components of the symmetric Cauchy stress tensor,𝕋, in a Walters-B fluid are given by [55–57]
(1)τij=-pδij+2∫0∞‍N(ξ)ξdξ∫-∞t‍e((t′-t)/ξ)∂xi∂x′m∂xj∂x′re(1)mrdt′,
where τij is the stress tensor, p represents the isotropic pressure, N(ξ)dξ is the total viscosity of the Maxwell elements whose relaxation time lies between ξ and dξ, (x′) represents the coordinates of a fluid element at time t′, (x) designates the coordinates of the same fluid element at t(>t′), δij is the metric tensor, and e(1)ij indicates the rate of strain tensor. Neglecting ∫0∞‍ξnN(ξ)dξ, n≥2 for liquids with short relaxation times, (1) reduces to a model involving only one non-Newtonian parameter [12, 57]:(2a)τij=-pδij+2μeij-2Γ0δeijδt,(2b)δeijδt=∂eij∂t+uk∂eij∂xk-∂uj∂xkeik-∂ui∂xkekj,where Γ0=∫0∞‍ξN(ξ)dξ is the viscoelastic parameter exhibiting low relaxation times, μ=∫0∞‍N(ξ)dξ is the limiting dynamic viscosity at small rates of shear, δ/δt evaluates the upper convected time derivative or Oldroyd derivative of a tensor [58], and ui are the components of the velocity u=(u,v,w) in the xi direction (x,y,z) of a fixed coordinate system. Refer to Appendix A for the complete expressions of the components of the stress tensor.Since the model proposed by Dandapat and Samanta [52] is for weakly viscoelastic fluids, it yields the same evolution equation for the free surface which results by adopting (2a)-(2b) [59]. This is true since β0 [52], the coefficient of cross-viscosity [60] defined by Dandapat and Samanta [52], is ignored by assuming the fluid to be weakly elastic, which occurs exactly by neglecting ∫0∞‍ξnN(ξ)dξ,n≥2.Such an idealized model represents an approximation to first order in elasticity [12] and is a valid approximation of Walters-B liquid taking very short memory [61]. Barnes et al. [62] have introduced a general canonical form of constitutive equations for viscoelastic fluids. The rheological model given by (2a) and (2b) has been recovered as a special case of this general form and is classified as a second-order model according to their nomenclature. Also, the justifications proposed by Dandapat and Samanta [52], Gupta [63] and Dandapat and Gupta [64] qualify and validate (2a)-(2b).The equations governing the flow are those for an incompressible fluid such that the compressibility effects and viscous heat generation are neglected in the heat equation:(3)ρ∂u∂t+ρu·∇u=∇·𝕋+ρF,(momentumbalance),∇·u=0,(continuityequation),∂T∂t+u·∇T=χ∇2T,(energybalance),
where F=(gsinβ,0,-gcosβ) is the body force, g is the standard gravity constant, and T is the temperature. To complete the mathematical description of (3), appropriate boundary conditions are prescribed at z=l(x,y) (bottom topography) and z=f(x,y,t) (interface location). They are
(4)u=0,T=TWonz=l(x,y),(no slip condition, uniform temperature distribution)(5)n·𝕋·n=-σ∇·nonz=f(x,y,t),(balance of normal stress at the interfacial surface)(6)ti·𝕋·n=ti·∇σonz=f(x,y,t),(balance of shear stress at the interfacial surface)(7)-λc∇T·n=α(T-T∞)onz=f(x,y,t),(Newton’s law of cooling at the free surface)(8)w=ht+u·∇fonz=f(x,y,t),(kinematic boundary condition at the free surface).In (5)–(7) above, n=∇(z-f)/|∇(z-f)| is the unit outward normal at any point on the free surface, and tx and ty are the unit tangential vectors such that n·ti=0, i=x,y. The heat exchange coefficient between the liquid and the passive air is described by α.In order to nondimensionalize systems (3)–(8), the problem is related to the flow of a falling film on an infinitely extending rigid plane [28, 38, 39]. For a stationary flow of uniform thickness h0, the solution of the Navier-Stokes and Fourier (expressing energy balance) equations is a parallel flow dependent only upon the cross-stream coordinate such that viscosity balances gravity and the heat is propagated by pure conduction. The primary flow variables [65] are (neglecting the capillary boundary layers at the side walls [66])
(9)u-=ρgsinβμ(zh0-z22),v-=0,w-=0,p-=p∞+ρgcosβ(h0-z),T-=TW-α(TW-T∞)λc+αh0z.
The depth averaged velocity of the Nusselt flow is defined as u0=(1/h0)∫0h0‍u-(z)dz = (ρgsinβh02)/3μ, and it serves to scale the streamwise component of the velocity, for which the constant flux is given by q0=u0h0. The Nusselt thickness is considered as the characteristic length scale in the z-direction. The characteristic wavelength on the free surface is considered to be of the same order as the wavelength of the periodic undulation and serves as a scale in the x- and y-directions. Mass balance condition serves to choose the scales for the velocity components along the spanwise and transverse directions. The pressure is scaled using ρu02 and the viscoelastic parameter [13] by ρh02. Lastly, the temperature scale is chosen to allow the largest possible temperature difference [21, 33]. In short, the scales can be summarized as
(10)(x*,y*,z*,h*,l*,f*)=1h0(εx,εy,z,h,l,f),t*=εu0th0,p*=p-p∞ρu02,(u*,v*,εw*)=1u0(u,v,w),T*=T-T∞TW-T∞,Γ*=Γ0ρh02,ε=h0λ≪1.After nondimensionalizing, the asterisks are dropped from the scaled variables (cf. AppendixB for the complete set of nondimensional equations). The set of dimensionless parameters arising during the nondimensionalization procedure is Re=ρu0h0/μ (Reynolds number), Γ=Γ0/ρh02 (viscoelastic parameter), We=σ∞/(ρu02h0) (Weber number), Ca=γ(TW-T∞)/σ∞(measures the ratio of temperature dependent part of the surface tension over its basic value; note that it is defined sometimes as the Capillary number [18, 19]), Pr=ν/χ (Prandtl number), Pe=RePr (Péclet number), M=γ(TW-T∞)/μu0 (Marangoni number), and Bi=αh0/λc(Biot number). Referring to Joo et al. [18] and the references therein, the orders of magnitude assigned to the dimensionless parameters are Re≃𝒪(1), ε2We=W≃𝒪(1), Ca≃𝒪(ε2), Γ≃𝒪(1), Pr≃𝒪(1), Pe≃𝒪(1), M≃𝒪(1) and Bi≃𝒪(1). Further, it is assumed that Al=al/h0≃𝒪(ε)(al/λ≃𝒪(ε2),h0/λ=ε) or Al=al/h0≃𝒪(1)(al/λ≃𝒪(ε),h0/λ=ε), where ε is the small aspect ratio. Note that the assumption Al≃𝒪(ε) allows to construct a series solution analytically as mentioned later [26, 30]. However, when Al≃𝒪(1), the base solution is found through a numerical procedure [26, 30].The physical variablesu, v, w, p, and T are asymptotically expanded in powers of the small parameter ε and solved sequentially at each order [21, 67] using the symbolic computation toolbox in MATLAB to obtain the long-wave model (cf. Appendix B for the equations and solutions at 𝒪(1) and 𝒪(ε)):
(11)ht+3h2hx+ε{Mh22∇(11+Bih)-5Re8∂∂x(h4ht)-27Re40∂∂x(h6hx)111111+32ReΓ∂∂x(h2ht)+152ReΓ∂∂x(h4hx)111111+∇·[Mh22∇(11+Bih)W3Reh3∇∇2f-cotβh3∇f111111111111-Mh22∇(11+Bih)]}+𝒪(ε2)=0.
It is to be noted that the above long-wave model is obtained without estimating the time derivatives from the leading order solutions as traditionally done [35, 36]. For small Reynolds numbers, the long-wave model described by (11) captures the dynamics efficiently since the temperature and velocity fields are slaved to the kinematics of the free surface. The second term in (11) describes wave propagation and steepening effect. The remainder of the terms in (11) within the curly braces decides the growth or decay of waves of which the first two terms are due to the inertial contribution and correspond to the destabilizing effects of the mean shear flow, the third and the fourth terms measure jointly the viscoelastic effects, and the fifth term is responsible for the effects of mean surface-tension. While the sixth term accounts for the stabilizing effects of the hydrostatic pressure, the seventh term is due to thermocapillarity and represents the perturbation of the interface temperature induced by variations of h, when there is heat transfer to the gas phase.


## 3. Temporal Linear Stability and Results

It is assumed that the primary flow,h=H(y), is independent of x and t and depends only upon the spanwise coordinate y. The primary flow satisfies the following ordinary differential equation:
(12)ReW3H2(d4ldy4+d4Hdy4)+ReWHdHdy(d3ldy3+d3Hdy3)-cotβH2(d2ldy2+d2Hdy2)-3cotβHdHdy(dldy+dHdy)+MBiH2(1+BiH)2d2Hdy2+MBi(1+BiH)3(dHdy)2=0.ForAl≪1 (a weak undulating topography), a closed form expression for the thickness profile of the primary flow can be computed [26, 30, 37–39] by expanding H(y) = 1+AlH1(y)+𝒪(Al2) such that H1(y)=C1sin(2πy)+C2cos(2πy) is periodic and represents the first order correction to the primary flow. This assumption used in (12) yields C1=0 and C2=-(C21/(C21+C22)), where C21=(2π)4(ReW/3)+(2π)2cotβ and C22=-(MBi/2(1+Bi)2). Since ε is small, another alternative to construct a mathematical form for the primary flow would be to set H(y)=1+εH1(y)+𝒪(ε2). This ansatz, however, imposes a mathematical condition on the Reynolds number, angle of inclination, wall amplitude and the Weber number associated with the flow at the leading order. And at 𝒪(ε), it involves solving a quartic auxiliary equation associated with the differential equation for determining H1(y), which as a result yields complicated complementary functions depending upon the nature of roots. Therefore, this approach is avoided here. Having known the undisturbed thickness of the basic state, the time-independent properties of the fluid like velocity, temperature, and so forth can be evaluated at the free surface and averaged over a typical groove [33] by referring to the leading and first order solutions in Appendix B. This yields(13a)u~=∫01‍(u0+εu1)|z=l+hdy=32+34Al2C22,(13b)T~=∫01‍(T0+εT1)|z=l+hdy=0.The above equations give a first insight on the flow mechanism associated with the topography comprising parallel grooves. From (13a), it is inferred that the averaged time-independent stream velocity enhances (which would cause the averaged flow rate, q~a, to increase) as the groove amplitude increases, suggesting possible flow destabilization. Also, the temperature at the bottom topography is larger than the averaged temperature at the free surface (cf. (B.2) and (13b), and note that the temperature averaged on a groove at the bottom substrate has the value 1). It is known that when the interfacial flow is driven by local changes in interfacial tension due to heating of the substrate, the modification of the temperature at the free surface due to surface deformation triggers long-wave thermocapillary instability [68, 69], sometimes leading to unwetted regions [70]. Therefore, the difference in temperature affects the surface tension and causes the fluid to flow from hotter region to colder regions, thereby increasing the wave amplitude.It should be remarked that contrary to the observations in Gambaryan-Roisman and Stephan [35], Sadiq et al. [37], and Sadiq [38], (13a) proposes possible flow destabilization on a patterned topography as the averaged flow rate over a typical groove increases. However, at this juncture, it is really premature to speculate and conclude on this behavior of the time-independent velocity on a grooved topography, given the fact that (13a) and (13b) are obtained only for small Al. Adopting a numerical solution of H(y) (as discussed later), a detailed investigation on the equilibria would only reveal the true picture of the basic state getting either stabilized or destabilized. In this regard, a small harmonic disturbance (a function of y times a sinusoid, which is a function of x and t) is imposed on the primary flow thickness to investigate the temporal stability and to see whether the disturbance dies away in the course of time (stable flow), persists as a disturbance of similar magnitude (neutral stability) or continues to grow in time leading to a different laminar or turbulent flow (instability) in the form,
(14)h=H(y)+ζη(y)eikx+λt,
where λ=CR-ikCL is the temporal pulsation (eigenvalue) such that k, CR, and CL are real numbers and represent the wavenumber, the linear amplification, and the linear wave velocity, respectively (by fixing k to be real, the solution can oscillate but is not allowed to grow/decay in x). The absolute value of the disturbance amplitude function ζη(y) is small such that max(|η(y)|)=𝒪(1) and ζ≪1. Substituting (14) in (11), a linear ordinary differential equation of fourth order in η(y) is found:
(15)ε3ReWH3d4ηdy4+εReWH2dHdyd3ηdy3+ε[MBiH22(1+BiH)2-cotβH3-23ReWk2H3]d2ηdy2+ε[+ReWH2(d3Hdy3+d3ldy3)-ReWk2H2dHdy]2MBi(1+BiH)3HdHdy-3cotβH2(2dHdy+dldy)11111+ReWH2(d3Hdy3+d3ldy3)-ReWk2H2dHdy]dηdy+ε[-3MBi2H(1+BiH)4(dHdy)2]3εikH2+2740k2ReH6-152ReΓk2H4111111+cotβk2H3+13ReWk4H3111111-6cotβHdHdy(dHdy+dldy)111111-3cotβH2(d2Hdy2+d2ldy2)111111+2ReWHdHdy(d3Hdy3+d3ldy3)111111+ReWH2(d4Hdy4+d4ldy4)111111-MBik2H22(1+BiH)2+MBiH(1+BiH)2d2Hdy2111111-MBi2H2(1+BiH)3d2Hdy2+MBi(1+BiH)3(dHdy)2111111-3MBi2H(1+BiH)4(dHdy)2]η=λ(-1-3iεkReΓ2H2+5iεkRe8H4)η,
provided that |eλt|<(1/Bi|ζη|)(1+Bi|H|).For the flow on a planar wall, the terms contributing to the real and imaginary parts of the complex eigenvalueλ are found as (H=1, l=0, and η is a constant)(16a)CR=(εk2[(65)Re+3ReΓ-cotβ+MBi2(1+Bi)2-ReW3k2])×[1+ε2k2(3ReΓ2-5Re8)]-1,(16b)CL=3+ε[3ReΓ2-5Re8]CR.For ε≪1, the above result reduces to the asymptotic analysis result of the Orr-Sommerfeld system considered in Sadiq et al. [21] for a Newtonian fluid, however on a rigid impermeable plane. The condition for the surface wave instability (CR>0) from (16a) becomes
(17)Re>cotβ-(MBi/2(1+Bi)2)6/5+3Γ-k2(W/3),
which shows that the critical Reynolds number decreases when the viscoelastic effect increases and also when the thermocapillary effect increases, indicating flow destabilization. However, when W increases the critical Reynolds number increases implying stabilization of the flow. Also, note that ignoring surface tension and heating effects at the leading order for a Newtonian film, the well-known hydrodynamic instability condition Re>(5/6)cotβ is recovered, which shows that a vertically falling film is always unstable.Whenη is a function of y and the flow is on a planar wall, (15) becomes a fourth-order linear differential equation with constant coefficients, which can be straightforwardly solved. On a topographical substrate when Marangoni and surface-tension effects are neglected, (15) reduces to the form p(y)(d2η/dy2)+q(y)η=0. For an eigenvalue λ, the eigenfunction η could be explicitly solved using a series solution technique for the above second-order differential equation at ordinary points y≠(1/2π)cos-1(1/Al).Since the coefficients of (15) are periodic with period 1, according to Floquet theory, η can be expressed in the form [71] η(y)=ϑ(y)eiςy, where ϑ(y) has period 1. The parameter ς is taken real by neglecting the exponential growth of disturbances along the spanwise direction. Together with the periodic boundary conditions for ϑ(y), when the topography is uneven, the eigenvalue problem (15) can be solved numerically to determine CR(k,ς) and CL(k,ς). Also, note that (15) remains invariant under the transformation y→-y, suggesting that (15) is symmetric about the (x,z) plane (remark that (11) is actually symmetric about this plane), and hence, one can concentrate only in the region y∈[0,1/2] and demand dη/dy=0, d3η/dy3=0 at the end points. Another alternative would be the application of symmetry boundary conditions along y=0 and periodic boundary conditions at the end points of a typical groove. However, in the rest of the analysis, only periodic conditions are assumed.In order to solve the linear stability problem (15), a spectral collocation method is used by considering periodic boundary conditions. The interval [0,1] in which the bottom topography l=Alcos(2πy) is defined and periodic in y is mapped onto [-1,1] using the linear transformation Y=2y-1. The collocation points Yj in the interval [-1,1] are distributed according to the cosine law (where N is a positive even integer):
(18)Yj=cos(jπN),j=0,1,2,…,N-1,N.
The derivatives d/dY are replaced by a matrix D of dimension (0:N,0:N), whose elements are(19a)D(r,s)=cr(-1)r+scs(Yr-Ys),0≤r,s≤N,r≠s,(19b)D(r,r)=-Yr2(1-Yr2),1≤r≤N-1,(19c)D(0,0)=-D(N,N)=2N2+16,where Yr,Ys are distributed according to the cosine law (18). Further, c0=cN=2, cs=1, and 1≤s≤N-1. The vector obtained by multiplying η with the matrix D is the discrete approximation of dη/dY. Higher order derivatives are approximated and obtained by repeated multiplication by the matrix D. For the cosine law based on the odd values of N and for the theoretical foundation for approximating the derivatives through successive matrix multiplications, the readers are advised to refer to the standard literature on spectral methods [72, 73].Instead of choosing a particular fluid and analyzing its stability (which has only a limited scope), the parameters are changed by varying the components ofP=(β, Al, Re, W, Γ, M, and Bi) to analyze the equilibria. However, following Joo et al. [18] and Miladinova et al. [19], ε=0.2 is fixed. Also, while varying the parameters such as Re and Bi, it is ensured that these parameters do not take unrealistic values. For instance, the long-wave equation on a planar surface displays finite-time blow-up behavior when the Reynolds number exceeds a particular value [20, 74], and for liquid films in contact with gases in many realistic experimental situations, Biot number is usually small [20]. Further, as the rheological model used in the study is assumed to be realistic only for weakly elastic fluids, the resulting flow should be considered as a perturbation of a Newtonian viscous flow, and therefore, small values of Γ are chosen. Lastly, appreciable surface-tension effect is included bearing in mind the thinness of the film.In the numerical calculations performed, the maximal values ofCR correspond to ς=0, so that the most unstable modes are those for which ϑ(y) is periodic. Therefore, the dependence of CR on k alone is considered from this point onwards. Note that the numerically generated solution of (12) was considered in the study for plotting all the figures presented here. The reason for doing this is twofold: (i) the accuracy of H(y) theoretically constructed using the imposition Al≪1 can be checked and (ii) solutions at Al~𝒪(1) can be computed. This was accounted by integrating (12) numerically by approximating the derivatives using the spectral collocation technique and then solving the resulting nonlinear algebraic system using Newton’s method. However, an initial condition is required to initiate integration. For this purpose, the evaluated analytical expression of the primary flow was considered as the initial condition, and periodic conditions were chosen as boundary conditions. Figure 2 displays the free surface obtained through approximated analytical and numerical solutions of H(y). It is observed that the numerical solution of H(y) shows a largest discrepancy of 8% from the analytical solution. This reduces further either when the magnitude of the Weber number increases or when the topography is more inclined, leading to agreement between the numerical and analytical values of Z(y), as shown in Figure 2. Interestingly, D’Alessio et al. [33] also reported on the behavior of free surface approaching unity as the magnitude of surface tension increases (cf. Figures4and23 in their research paper). Finally, note that in order to ensure the accuracy of the numerical results yielded by the spectral method, the derivatives were approximated using a central difference scheme of second-order accuracy and solved (only the numerical results obtained through spectral method are however presented in this paper). The numerical results ensuing due to these two different procedures agreed with each other. For this purpose, N≥40 was considered in the spectral technique, and N≥200 was chosen to implement the finite difference scheme.Free surface profiles corresponding to the primary flow on an uneven substrate withAl=0.5, Re=5, M=5, and Bi=0.1: (a) cotβ=2, W=1; (b) cotβ=2, W=10, and (c) cotβ=0, W=1.
(a)
(b)
(c)Stability characteristics of Newtonian and non-Newtonian isothermal layers are examined first. It is known that the flow of an isothermal liquid film with a smooth free surface down a vertical plane is unstable at all Reynolds numbers [75]. Also, since the amplitude of the long waves on a vertically falling film is larger than that along an inclined substrate case [37, 38], in the rest of the investigation, β=π/2 is set.It is ascertained from the numerical behavior of the linear amplification curves that the inertial force increases the instability threshold for Newtonian films falling along parallel longitudinal grooves (Figure3). For weakly undulating topography at small Reynolds numbers, the grooves aligned parallel to the flow direction decrease the growth rate leading to waves with small amplitude, and reduces the range of unstable wavenumbers in relation to the flow on a flat substrate. As the Reynolds number increases, there exists a certain wavenumber below which the growth rate for the flow along grooves is larger than that along a planar substrate, and such a trend changes beyond this wavenumber (Figure 3(b) (i) and Figure 3(c) (i)). Take a note on the maximum amplitude profiles presented in Sadiq et al. [37] and Sadiq [38], where Hmax|Al>0>Hmax|Al=0 for a short period of time before this behavior changes. The probable explanation to this behavior can be gathered from (13a), which suggests increase in average velocity at the free surface on a typical groove. Therefore, when a small disturbance is imposed on the primary flow, the grooves tend to increase the growth rate initially. When the bottom undulations increase, the numerical investigation exhibits interesting behavior on the films stability. There are instances at which the maximal growth rate corresponding to a large Weber number is larger than that related to a small Weber number; however, the range of unstable wavenumbers remains smaller than the case corresponding to small Weber numbers (Figure 3(b) (ii)). For a fixed moderate groove amplitude, when the Weber number increases, the flow destabilizes at small Reynolds numbers (Figure 3(a) (ii) and Figure 3(b) (iii)) and stabilizes when the magnitude of the inertial force increases (Figure 3(c) (iii)). Pointing out on the observations reported through direct [27, 28, 33] and indirect problems [76] for the flow across wavy undulations on an inclined topography (where the critical Reynolds number decreased on an undulating topography when the effect of surface tension is appreciable), here, at fixed Re and Al, the variation in We shows either destabilizing or stabilizing tendency on the flow.Variation of growth rate with wavenumber for a vertically falling Newtonian isothermal film: (a)Re=0.5, (b) Re=1 and (c) Re=5; (1)W=1, (2)W=5, and (3)W=10. In a(i)–c(i), solid lines represent Al=0, and broken lines are the results for Al=0.4, and (ii) and (iii) are for Al=0.7 and Al=0.9, respectively.
(a)
(b)
(c)When the non-Newtonian effects are included, the infinitesimal disturbance to the primary flow leads the viscoelastic film to destabilize in relation to the Newtonian film. The sketch in Figure4 displays interesting trends associated with the growth rate curves of isothermal Newtonian and viscoelastic layers on a vertical substrate, when there is a pronounced inertial effect. At Re=5, the numerical results of the stability analysis predict decrement in the value of growth rate when the surface-tension effect increases (compare the curves in each column in Figure 4). Analogous to the Newtonian flow, the growth rate curves corresponding to the viscoelastic film on a grooved topography show increasing and decreasing behavior before and after a certain wavenumber in connection to the flow on a flat substrate. However, at large Weber numbers, when the viscoelastic parameter increases, the growth rate and the window of unstable wavenumbers for the flow on a topographical substrate are larger than the case when Al=0. For example, referring to Figure 4(c) (ii) and Figure 4(c) (iii), it is seen that Al=0.4 yields larger growth rate values than Al=0, and the range of unstable wavenumbers in Figure 4(c) (ii) corresponding to Al=0 (unstable k∈[0,0.933]) is smaller than that for Al=0.4 (k∈[0,0.9445]). For Figure 4(c) (iii), unstable k for Al=0 lies in the interval [0,0.664], and that for Al=0.4 lies in [0,0.6747]. But when the amplitude of the undulating groove becomes larger, the growth rate curves and the range of unstable wavenumbers start decreasing (cf. curve 4 in Figure 4(c) (ii) and Figure 4(c) (iii)). See also Figure 11 presented later displaying the above characteristics on a heated uneven substrate. Figure 5 shows the dispersion curves for linear growth rate of a viscoelastic film at small Reynolds numbers. Just as in the case of a Newtonian film, at fixed moderate groove amplitudes, the topographical substrate destabilizes the small Reynolds number viscoelastic flow as the magnitude of the surface-tension effect increases (Figure 5 (ii), (iii)). As a note from the numerical study carried, it should be mentioned that when a value as high as Γ=1 was considered (such a value may be beyond the realistic range of weakly elastic fluids), the longitudinal grooves destabilized the flow completely.Dispersion of linear waves for isothermal Newtonian and viscoelastic films, temporal growth rate atRe=5, (a) Γ=0, (b) Γ=0.02 and (c) Γ=0.06; (1)Al=0, (2)Al=0.4, (3)Al=0.7, and (4)Al=0.9; (i) W=1, (ii) W=5, and (iii) W=10.
(a)
(b)
(c)Variation of growth rate with wavenumber for a vertically falling viscoelastic isothermal film withΓ=0.06 and Re=0.5, (1)W=1, (2)W=5, and (3)W=10. In (i), solid lines represent Al=0, and broken lines are the results for Al=0.4, and (ii) and (iii) are for Al=0.7 and Al=0.9, respectively.
(a)
(b)
(c)To substantiate the information availed from Figures3, 4, and 5, neutral stability curves are presented. While the neutral stability curves in Figure 6 clearly display the destabilizing nature of the periodic bottom undulations at small Reynolds and large Weber numbers, Figure 7 distinctly shows the instability at small Al when the viscoelastic effects increase at Re=5, as already discussed in relation to Figure 4.Stability chart ink-Al plane for a vertically falling isothermal film: (a) Γ=0, (b) Γ=0.06; (i) Re=0.5, (ii) Re=1, and (1)W=1, (2)W=5, and (3)W=10 10; S-Stable region, US-Unstable region.
(a)
(b)Neutral stability curves of an isothermal film atRe=5 and β=90∘: (a) W=5, (b) W=10, and (1)Γ=0, (2)Γ=0.02, and (3)Γ=0.06.
(a)
(b)In Figure8, the phase velocity curves varying with k are displayed. The phase speed value is 3 for the flow on a planar substrate initially (this agrees with the theoretical prediction suggested by (16b) at ε≪1). However, the nonlinear nature of the phase velocity curves can be attributed due to the gross effect of all the dimensionless parameters involved in the problem without truncating or approximating them at any level; for example, see (16b). The linear wave velocity curves first decrease up to a certain wavenumber and then increase beyond it. It is seen that the phase velocity decreases as the viscoelastic parameter increases. This fact can be confirmed by referring to the flow over a flat substrate (cf. (16b)), where it is observed that when CR>0, CL would increase as a function of Γ only when Γ>5/12. Otherwise, despite CR>0, CL can decrease when Γ increases provided that 0≤Γ≤5/12. Only numerical calculations reveal the influence of an uneven geometry on the wave velocity. Computations show that the phase velocity is larger for the case of a flow on a grooved topography than on a planar substrate.Figure 8
Phase velocity curves:(1)P(π/2,0,5,5,0,0,0), (2)P(π/2,0.4,5,5,0,0,0), (3)P(π/2,0,5,5,0.02,0,0), (4)P(π/2,0.4,5,5,0.02,0,0), (5)P(π/2,0,5,5,0.06,0,0), and (6)P(π/2,0.04,5,5,0.06,0,0).For the case of a topographical substrate subjected to uniform heating, D’Alessio et al. [33] reported on very interesting phenomena associated with Marangoni effect. Their study witnessed that thermocapillarity can either stabilize or destabilize the flow depending upon the strength of the Marangoni number and can also lead to a reversal in flow stability. Figure 9 shows the response of thermocapillary effect on Newtonian and viscoelastic films over a flat and uneven topography. At low Reynolds number, the amplitude of the grooves decides the instability mechanism depending upon the strength of heating. For small groove amplitudes and for a fixed Weber number, increase in Marangoni number leads to film destabilization. Apart from the already observed behavior of flow destabilization at moderate reasonable Al with increase in Weber number, it is seen that increase in Marangoni effect can promote stability at large Weber numbers (cf. Figure 9(c)). Such an effect is also observed when the Reynolds number increases (figures are not included here in order to avoid repetition). Observe that in the recent research report of D’Alessio et al. [33], they have identified on a semiinclined plane (cf. Figure 8 in their paper) that for a given groove measure with increase in Marangoni number, the critical Reynolds number increases leading to film stabilization, provided the surface-tension effect is large. When the viscoelastic effects are included and when the Reynolds number associated with the flow is very small, the film tends to destabilize in the absence of Marangoni effect. However, at Al=0.7 and for pronounced surface-tension effects, increase in heating corresponds to flow stabilization (in Figure 9(d), compare curves 3 with 4 and 5 with 6). The influence of Marangoni stress on the stability threshold at low Reynolds number is presented by means of Figure 10, which shows improved stability at large Weber numbers at certain regions of Al. However, it is also inferred from the figure that as Al increases, the trend changes, and once again Marangoni stresses start destabilizing the flow. In conjunction with the results displayed in Figure 4(c) (ii), Figure 11 reveals film destabilization at Re=5 when the viscoelastic effects increase on a grooved topography (cf. curves for Al=0 and 0.4 in Figure 11). The effect of thermocapillarity in this case is to increase the growth rate. However, at moderate Al=0.9 increase in heating stabilizes the flow, similar to the case discussed in Figure 9(d).The growth rate of the amplitudes atRe=0.5 and Bi=0.1 on a vertically held substrate. The curves in (a), (b), and (c) correspond to Newtonian films. In (d), the odd numbers represent M=0, and the even numbers are for M=10. Also, the solid lines in (d) are for Γ=0, the broken lines are for Γ=0.02, and the broken-dotted lines are for Γ=0.06. Note that although Bi=0.1 is taken, in case of M=0, it corresponds to an isothermal problem, and the effect of Bi is ignored.
(a)
(b)
(c)
(d)Figure 10
Influence of Marangoni stress on the neutral stability curves atRe=0.5 and Bi=0.1 on a vertically held substrate: (1)W=5 and (2)W=10. Broken lines are for M=0, and the solid lines represent M=10.Figure 11
Growth rate as a function of wavenumber atRe=5 and Bi=0.1.Figure12 presents linear amplification curves at Re=5 for different values of Biot number for both Newtonian and non-Newtonian liquids. When the rate of heat transfer from the liquid to the ambient gas phase increases, the flow tends to destabilize (Figures 12(a) and 12(b)). For a non-Newtonian film on a planar substrate with increase in the magnitude of Bi, the growth rate curves dominate the Newtonian case causing more flow destabilization (Figure 12(a)). On a substrate with topography, the groove amplitude plays a crucial role in deciding films stability. Note that increase in the value of viscoelastic parameter for a flow on a topography for which Al=0.4 causes flow destabilization compared to the planar case at small values of Bi as inferred from the curves 4 in Figures 12(a) and 12(b) for which the range of unstable wavenumbers, respectively is [0, 0.984] and [0, 0.972], despite the respective maximal growth rates being CRmax=0.27 and CRmax=0.303. However, when the heat transfer rate increases, the growth rate curves over a topography with small amplitude increase up to a certain wavenumber in relation to the growth rate curve over a flat substrate, and this trend changes beyond such a wavenumber causing reversal in stability (cf. curves 6(a) and 6(b) in Figure 12(d), where these curves correspond to curve 6 from Figures 12(a) and 12(b), resp.). At moderate values of groove amplitude, it was already observed that the Marangoni effect can stabilize the flow mechanism. This stabilizing mechanism enhances further when the heat transfer rate increases (Figure 12(c)). When the Reynolds number corresponding to the flow is small and when W and M are large, increase in Biot number stabilizes the flow configuration as seen in Figure 13 (notice that the stabilization process corresponding to increased thermocapillary effect is also seen).Linear amplification rate:(1)P(π/2,Al,5,5,0,0,0), (2)P(π/2,Al,5,5,0.06,0,0), (3)P(π/2,Al,5,5,0,5,0.1), (4)P(π/2,Al,5,5,0.06,5,0.1), (5)P(π/2,Al,5,5,0,5,0.5) and (6)P(π/2,Al,5,5,0.06,5,0.5); (a) Al=0, (b) Al=0.4, and (c) Al=0.9. In (d), the curves correspond to (2) and (6) from (a) and (b).
(a)
(b)
(c)
(d)Figure 13
Curves corresponding to the growth rate on a vertical topography withAl=0.7.The present investigation, although valid only at the region where the viscous forces dominate inertia due to long-wave modeling, has revealed interesting results associated with the linear stability. The results showed the dramatic impact of the bottom topography on the films stability in conjunction with the force of surface-tension and Marangoni stresses. The analytical expansion of the streamwise velocity at the free surface averaged over a typical bottom undulation does reveal only its increasing tendency contributing towards flow destabilization, but the numerical investigation on the films stability exposes the additional effects involved in the problem. A close-up theoretical investigation can be performed by referring to Figure2 to get an insight into the physical mechanism involved. In this figure, it is observed that the free surface flattens up at z=1 for a large value of W and a small value of bottom amplitude. Using this fact and referring to Appendix B, the leading order solutions can be obtained. For the streamwise flow, amongst others, the equations
(20)u0zz+3=0,u0=0onz=l(x,y),u0z=0onz=1
lead to
(21)u0(y,z)=3(z-l)+32(l2-z2).
An integration of (21) yields the base flow rate as
(22)q~=∫l1‍u0(y,z)dz=(1-l)3
and the averaged flow rate as
(23)q~a=∫01‍q~dy=1+3Al22,
indicating that for large values of surface tension force the averaged value of the flow rate increases leading to flow destabilization. The undisturbed thickness being a function of the spanwise variable, when perturbed along the x-direction, leads to an eigenvalue problem comprising terms arising as a consequence of streamwise perturbation. The net effect of these terms with the aid of numerical approach shows that both flow stabilization and destabilization are possible in various parametric regimes.


## 4. Concluding Remarks and Objective

The linear stability of a thin Walters-B fluid due to the linear variation of the surface-tension along a liquid interface induced by temperature gradients along the liquid layer was examined on a topography comprising sinusoidal longitudinal grooves. The nontrivial eigenvalue problem was numerically solved using a Chebyshev spectral collocation technique and the results were compared with those occurring from a central difference scheme of second order accuracy to ensure numerical correctness.Numerical investigation of the eigenvalue problem corresponding to a falling film on a vertically held substrate revealed interesting results. In principle, when the viscous effects dominate the inertial forces, the flow tends to stabilize the system on a longitudinal groove topography, provided the surface-tension effects are weak. However, when the magnitude of the surface-tension force increases, interesting trends were exhibited by the stability curves. Some of the interesting key observations in the study can be enlisted as follows: (i) increased surface-tension effects contribute towards flow destabilization at low-Reynolds number as the amplitude of the bottom undulations increase (ii) Marangoni stresses can lead the system towards flow stabilization for certain increasing range overAl, provided the strength of the Weber number is appreciable, (iii) viscoelastic effects destabilize the system compared to a Newtonian flow, but the structured topography tends to stabilize the flow in relation to the flow on a planar substrate for weak viscoelastic effects, (iv) increase in viscoelasticity promotes instability for small groove amplitudes when the Weber and Reynolds number magnitude increases and (v) increase in heat transport rate at large Weber numbers causes flow stabilization at low-Reynolds numbers.Calculations were performed through numerical approach to analyze the effect of perturbation on the undisturbed state on a longitudinal grooved topography. Theoretically, it was also found that the average flow rate increases when the surface-tension force increases at small groove amplitudes. Although the above facts are findings based on a numerical linear stability study, such a behavior is yet to be confirmed through laboratory investigations. However, the latest research article by D’Alessio et al. [33] reported strikingly on similar characteristics for the flow of a Newtonian film but across the undulations over an inclined topography.It is not only with a fundamental research point of interest but also due to benefits gained through peer research that the thin film community intensively and actively strives to shed light on understanding intricate and complex mechanisms involved in falling film problems and fill up unknown gaps. It should be remarked that Mazouchi and Homsy [77] compared the solution to two dimensional Stokes flow to their earlier thin-film based solution. This comparison indicated that the thin-film hypothesis remains valid even in the presence of steep topographic variation. In this regard, theoretical understanding and modeling of flows and deformation in Newtonian/non-Newtonian fluids on small scales on a patterned topography assuming the film to be thin would help in the design and fabrication of small scale structures over the topography to regulate heat and mass transfer in soft matter with smart and advanced functionalities to suppress instabilities and in preventing rupture mechanism. In addition, such studies may help in understanding the unwanted influence of surface irregularities, which result from a particular stage in a manufacturing process thereby affecting the flow mechanism.Walters-B fluid finds its usage in medicinal diagnosis and clinical applications. For instance, Nadeem et al. [78] have mentioned that viscoelastic fluid can easily flow in intestine because elastic materials strain instantaneously when stretched and just as quickly return to their original state once the stress is removed and that the Walters-B model is quite close to our real life system. The present study may serve as a part of the basis to study the peristaltic flow of Walters-B model in asymmetrical channels or tubes with geometries inside. Also, Walters-B type finds relevance and importance [60] in geophysical fluid dynamics, chemical technology, and petroleum industries. An implication of the study could be possibly in the product quality and patterning of the substrates, which require multilayer coating where one of the purposes of the bottom layer is to planarize the topography, that is, to smooth out the uneven topography. And a nonuniform coating leads to low quality products or to manufacturing failures. Pavlidis et al. [45] reported on the steady viscoelastic film flow over 2D topography with the prime objective to analyze the effect of the viscoelastic properties of polymeric solutions in spreading over and planarizing various topographical features, especially during spin-coating. Such studies are interesting and could be extended for topographies with longitudinal geometries for second-order fluids with or without invoking lubrication approximation.Mathematical models obtained by combining a systematic gradient expansion with weighted residual techniques using polynomials as test functions may possibly shed more light on the stability thresholds beyond the low-Reynolds number regime [20, 79]. Such a study has been accomplished by Amatousse et al. [80] to study traveling waves on falling Walters-B viscoelastic fluid. As a future objective, a Galerkin strategy will be employed to obtain a system of equations coupling the film thickness and flow rates at various orders to investigate the stability [20, 38, 79, 80]. Also, it would be interesting to investigate the combined effect of nonuniform heating and the substrate topography on the equilibria of a film, for the case when the flow rate is high and when the topography is inundated. Of particular interest would also be to analyze the phenomena in case of a partially filled topography occurring due to low flow rate, where the peaks of the topography are exposed and remain dry. Another direction of study could be to use an Oldroyd-B model and investigate the stability if the nonlinearities inherent in more realistic constitutive equations become important [47]. Such studies will be addressed in future investigations.



---

*Source: 751909-2013-11-20.xml*