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

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.

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×10³Kg/m³, limiting viscosity 𝜇=0.79Ns/m², kinematic viscosity 𝜈=806×10⁻⁶Nsm/Kg, surface tension 𝜎=40×10⁻³N/m, and viscoelastic coefficient Γ₀=0.04Ns²/m². 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 of 1<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 film across 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

Figure 1 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 𝑥-direction, and the geometry extends infinitely along the 𝑥- and 𝑦-directions. The bottom topography, 𝑙(𝑥,𝑦)=𝑎_𝑙cos(2𝜋𝑦/𝜆) undulates weakly along the 𝑦-axis such that 𝑓(𝑥,𝑦,𝑡)=𝑙(𝑥,𝑦)+ℎ(𝑥,𝑦,𝑡) measures the total thickness. The parameters 𝑎_𝑙 and 𝜆 represent the amplitude (the peak deviation of 𝑙 from its center position) and the wavelength of a typical periodic groove. The air above the liquid is passive such that its temperature, 𝑇_∞, and pressure, 𝑝_∞, are constants sufficiently far from the interface. A uniform temperature distribution 𝑇_𝑊 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 (𝜆_𝑐) and thermal diffusivity (𝜒) are assumed constants [1, 20, 48]. The surface tension decreases linearly as a function of temperature and is given by 𝜎 = 𝜎_∞−𝛾(𝑇−𝑇_∞), where 𝜎_∞ is the mean surface tension at gas temperature and 𝛾=−𝑑𝜎/𝑑𝑇 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] 𝜏_𝑖𝑗=−𝑝𝛿_𝑖𝑗+2^∞₀𝑁(𝜉)𝜉𝑑𝜉^𝑡_−∞𝑒^{((𝑡′−𝑡)/𝜉)}𝜕𝑥^𝑖𝜕𝑥^′𝑚𝜕𝑥^𝑗𝜕𝑥^′𝑟𝑒^(1)𝑚𝑟𝑑𝑡^′, (1) where 𝜏_𝑖𝑗 is the stress tensor, 𝑝 represents the isotropic pressure, 𝑁(𝜉)𝑑𝜉 is the total viscosity of the Maxwell elements whose relaxation time lies between 𝜉 and 𝑑𝜉, (𝑥^′) represents the coordinates of a fluid element at time 𝑡^′, (𝑥) designates the coordinates of the same fluid element at 𝑡(>𝑡^′), 𝛿_𝑖𝑗 is the metric tensor, and 𝑒^(1)𝑖𝑗 indicates the rate of strain tensor. Neglecting ∫^∞₀𝜉^𝑛𝑁(𝜉)𝑑𝜉, 𝑛≥2 for liquids with short relaxation times, (1) reduces to a model involving only one non-Newtonian parameter [12, 57]:𝜏_𝑖𝑗=−𝑝𝛿_𝑖𝑗+2𝜇𝑒_𝑖𝑗−2Γ₀𝛿𝑒_𝑖𝑗,𝛿𝑡𝛿𝑒_𝑖𝑗=𝛿𝑡𝜕𝑒_𝑖𝑗𝜕𝑡+𝑢_𝑘𝜕𝑒_𝑖𝑗𝜕𝑥_𝑘−𝜕𝑢_𝑗𝜕𝑥_𝑘𝑒_𝑖𝑘−𝜕𝑢_𝑖𝜕𝑥_𝑘𝑒_𝑘𝑗, (2a)(2b)where Γ₀=∫^∞₀𝜉𝑁(𝜉)𝑑𝜉 is the viscoelastic parameter exhibiting low relaxation times, 𝜇=∫^∞₀𝑁(𝜉)𝑑𝜉 is the limiting dynamic viscosity at small rates of shear, 𝛿/𝛿𝑡 evaluates the upper convected time derivative or Oldroyd derivative of a tensor [58], and 𝑢_𝑖 are the components of the velocity 𝐮=(𝑢,𝑣,𝑤) in the 𝑥_𝑖 direction (𝑥,𝑦,𝑧) 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 𝛽₀ [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 ∫^∞₀𝜉^𝑛𝑁(𝜉)𝑑𝜉,𝑛≥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: 𝜌𝜕𝐮(),,𝜕𝑡+𝜌𝐮⋅∇𝐮=∇⋅𝕋+𝜌𝐅,momentumbalance∇⋅𝐮=0,continuityequation𝜕𝑇𝜕𝑡+𝐮⋅∇𝑇=𝜒∇²𝑇,energybalance, (3) where 𝐅=(𝑔sin𝛽,0,−𝑔cos𝛽) is the body force, 𝑔 is the standard gravity constant, and 𝑇 is the temperature. To complete the mathematical description of (3), appropriate boundary conditions are prescribed at 𝑧=𝑙(𝑥,𝑦) (bottom topography) and 𝑧=𝑓(𝑥,𝑦,𝑡) (interface location). They are 𝐮=0,𝑇=𝑇_𝑊,,()𝐭on𝑧=𝑙𝑥,𝑦noslipcondition,uniformtemperaturedistribution𝐧⋅𝕋⋅𝐧=−𝜎∇⋅𝐧on𝑧=𝑓𝑥,𝑦,𝑡balanceofnormalstressattheinterfacialsurface_𝑖⋅𝕋⋅𝐧=𝐭_𝑖,()⋅∇𝜎on𝑧=𝑓𝑥,𝑦,𝑡balanceofshearstressattheinterfacialsurface−𝜆_𝑐∇𝑇⋅𝐧=𝛼𝑇−𝑇_∞on𝑧=𝑓𝑥,𝑦,𝑡,Newton'slawofcoolingatthefreesurface𝑤=ℎ_𝑡,.+𝐮⋅∇𝑓on𝑧=𝑓𝑥,𝑦,𝑡kinematicboundaryconditionatthefreesurface (4)(5)(6)(7)(8)

In (5)–(7) above, 𝐧=∇(𝑧−𝑓)/|∇(𝑧−𝑓)| is the unit outward normal at any point on the free surface, and 𝐭_𝑥 and 𝐭_𝑦 are the unit tangential vectors such that 𝐧⋅𝐭_𝑖=0, 𝑖=𝑥,𝑦. 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 ℎ₀, 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]) 𝑢=𝜌𝑔sin𝛽𝜇𝑧ℎ₀−𝑧²2,𝑣=0,𝑤=0,𝑝=𝑝_∞+𝜌𝑔cos𝛽ℎ₀−𝑧,𝑇=𝑇_𝑊−𝛼𝑇_𝑊−𝑇_∞𝜆_𝑐+𝛼ℎ₀𝑧. (9) The depth averaged velocity of the Nusselt flow is defined as 𝑢₀=(1/ℎ₀)∫^ℎ0₀𝑢(𝑧)𝑑𝑧 = (𝜌𝑔sin𝛽ℎ²₀)/3𝜇, and it serves to scale the streamwise component of the velocity, for which the constant flux is given by 𝑞₀=𝑢₀ℎ₀. The Nusselt thickness is considered as the characteristic length scale in the 𝑧-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 𝑥- and 𝑦-directions. Mass balance condition serves to choose the scales for the velocity components along the spanwise and transverse directions. The pressure is scaled using 𝜌𝑢²₀ and the viscoelastic parameter [13] by 𝜌ℎ²₀. Lastly, the temperature scale is chosen to allow the largest possible temperature difference [21, 33]. In short, the scales can be summarized as 𝑥^∗,𝑦^∗,𝑧^∗,ℎ^∗,𝑙^∗,𝑓^∗=1ℎ₀,𝑡𝜀𝑥,𝜀𝑦,𝑧,ℎ,𝑙,𝑓_∗=𝜀𝑢₀𝑡ℎ₀,𝑝^∗=𝑝−𝑝_∞𝜌𝑢²₀,𝑢^∗,𝑣^∗,𝜀𝑤^∗=1𝑢₀(𝑢,𝑣,𝑤),𝑇^∗=𝑇−𝑇_∞𝑇_𝑊−𝑇_∞,Γ^∗=Γ₀𝜌ℎ²₀ℎ,𝜀=⁰𝜆≪1. (10)

After nondimensionalizing, the asterisks are dropped from the scaled variables (cf. Appendix B for the complete set of nondimensional equations). The set of dimensionless parameters arising during the nondimensionalization procedure is Re=𝜌𝑢₀ℎ₀/𝜇 (Reynolds number), Γ=Γ₀/𝜌ℎ²₀ (viscoelastic parameter), We=𝜎_∞/(𝜌𝑢²₀ℎ₀) (Weber number), Ca=𝛾(𝑇_𝑊−𝑇_∞)/𝜎_∞(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), 𝑀=𝛾(𝑇_𝑊−𝑇_∞)/𝜇𝑢₀ (Marangoni number), and Bi=𝛼ℎ₀/𝜆_𝑐(Biot number). Referring to Joo et al. [18] and the references therein, the orders of magnitude assigned to the dimensionless parameters are Re≃𝒪(1), 𝜀²We=𝑊≃𝒪(1), Ca≃𝒪(𝜀²), Γ≃𝒪(1), Pr≃𝒪(1), Pe≃𝒪(1), 𝑀≃𝒪(1) and Bi≃𝒪(1). Further, it is assumed that 𝐴_𝑙=𝑎_𝑙/ℎ₀≃𝒪(𝜀)(𝑎_𝑙/𝜆≃𝒪(𝜀²),ℎ₀/𝜆=𝜀) or 𝐴_𝑙=𝑎_𝑙/ℎ₀≃𝒪(1)(𝑎_𝑙/𝜆≃𝒪(𝜀),ℎ₀/𝜆=𝜀), where 𝜀 is the small aspect ratio. Note that the assumption 𝐴_𝑙≃𝒪(𝜀) allows to construct a series solution analytically as mentioned later [26, 30]. However, when 𝐴_𝑙≃𝒪(1), the base solution is found through a numerical procedure [26, 30].

The physical variables 𝑢, 𝑣, 𝑤, 𝑝, and 𝑇 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 𝒪(𝜀)): ℎ_𝑡+3ℎ²ℎ_𝑥+𝜀−5Re8𝜕ℎ𝜕𝑥⁴ℎ_𝑡−27Re𝜕40ℎ𝜕𝑥⁶ℎ_𝑥+32𝜕ReΓℎ𝜕𝑥²ℎ_𝑡+152𝜕ReΓℎ𝜕𝑥⁴ℎ_𝑥𝑊+∇⋅3Reℎ³∇∇²𝑓−cot𝛽ℎ³−∇𝑓𝑀ℎ²2∇1𝜀1+Biℎ+𝒪²=0. (11) 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 ℎ, when there is heat transfer to the gas phase.

3. Temporal Linear Stability and Results

It is assumed that the primary flow, ℎ=𝐻(𝑦), is independent of 𝑥 and 𝑡 and depends only upon the spanwise coordinate 𝑦. The primary flow satisfies the following ordinary differential equation: Re𝑊3𝐻²𝑑⁴𝑙𝑑𝑦⁴+𝑑⁴𝐻𝑑𝑦⁴+Re𝑊𝐻𝑑𝐻𝑑𝑑𝑦³𝑙𝑑𝑦³+𝑑³𝐻𝑑𝑦³−cot𝛽𝐻²𝑑²𝑙𝑑𝑦²+𝑑²𝐻𝑑𝑦²−3cot𝛽𝐻𝑑𝐻𝑑𝑦𝑑𝑙+𝑑𝑦𝑑𝐻+𝑑𝑦𝑀Bi𝐻2(1+Bi𝐻)²𝑑²𝐻𝑑𝑦²+𝑀Bi(1+Bi𝐻)³𝑑𝐻𝑑𝑦²=0. (12)

For 𝐴_𝑙≪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 𝐻(𝑦) = 1+𝐴_𝑙𝐻₁(𝑦)+𝒪(𝐴²_𝑙) such that 𝐻₁(𝑦)=𝐶₁sin(2𝜋𝑦)+𝐶₂cos(2𝜋𝑦) is periodic and represents the first order correction to the primary flow. This assumption used in (12) yields 𝐶₁=0 and 𝐶₂=−(𝐶₂₁/(𝐶₂₁+𝐶₂₂)), where 𝐶₂₁=(2𝜋)⁴(Re𝑊/3)+(2𝜋)²cot𝛽 and 𝐶₂₂=−(𝑀Bi/2(1+Bi)²). Since 𝜀 is small, another alternative to construct a mathematical form for the primary flow would be to set 𝐻(𝑦)=1+𝜀𝐻₁(𝑦)+𝒪(𝜀²). 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 𝐻₁(𝑦), 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𝑢=¹₀𝑢₀+𝜀𝑢_{1𝑧=𝑙+ℎ}3𝑑𝑦=2+34𝐴²_𝑙𝐶²₂,𝑇=¹₀𝑇₀+𝜀𝑇_{1𝑧=𝑙+ℎ}𝑑𝑦=0. (13a)(13b)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, 𝑞_𝑎, 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 𝐴_𝑙. Adopting a numerical solution of 𝐻(𝑦) (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 𝑦 times a sinusoid, which is a function of 𝑥 and 𝑡) 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) where 𝜆=𝐶_𝑅−𝑖𝑘𝐶_𝐿 is the temporal pulsation (eigenvalue) such that 𝑘, 𝐶_𝑅, and 𝐶_𝐿 are real numbers and represent the wavenumber, the linear amplification, and the linear wave velocity, respectively (by fixing 𝑘 to be real, the solution can oscillate but is not allowed to grow/decay in 𝑥). The absolute value of the disturbance amplitude function 𝜁𝜂(𝑦) is small such that max(|𝜂(𝑦)|)=𝒪(1) and 𝜁≪1. Substituting (14) in (11), a linear ordinary differential equation of fourth order in 𝜂(𝑦) is found: 𝜀3Re𝑊𝐻³𝑑⁴𝜂𝑑𝑦⁴+𝜀Re𝑊𝐻²𝑑𝐻𝑑𝑑𝑦³𝜂𝑑𝑦³+𝜀𝑀Bi𝐻²2(1+Bi𝐻)²−cot𝛽𝐻³−23Re𝑊𝑘²𝐻³𝑑²𝜂𝑑𝑦²+𝜀2𝑀Bi(1+Bi𝐻)³𝐻𝑑𝐻𝑑𝑦−3cot𝛽𝐻²2𝑑𝐻+𝑑𝑦𝑑𝑙𝑑𝑦+Re𝑊𝐻²𝑑³𝐻𝑑𝑦³+𝑑³𝑙𝑑𝑦³−Re𝑊𝑘²𝐻²𝑑𝐻𝑑𝑦𝑑𝜂3𝑑𝑦+𝜀𝜀𝑖𝑘𝐻²+27𝑘40²Re𝐻⁶−152ReΓ𝑘²𝐻⁴+cot𝛽𝑘²𝐻³+13Re𝑊𝑘⁴𝐻³−6cot𝛽𝐻𝑑𝐻𝑑𝑦𝑑𝐻+𝑑𝑦𝑑𝑙𝑑𝑦−3cot𝛽𝐻²𝑑²𝐻𝑑𝑦²+𝑑²𝑙𝑑𝑦²+2Re𝑊𝐻𝑑𝐻𝑑𝑑𝑦³𝐻𝑑𝑦³+𝑑³𝑙𝑑𝑦³+Re𝑊𝐻²𝑑⁴𝐻𝑑𝑦⁴+𝑑⁴𝑙𝑑𝑦⁴−𝑀Bi𝑘²𝐻²2(1+Bi𝐻)²+𝑀Bi𝐻(1+Bi𝐻)²𝑑²𝐻𝑑𝑦²−𝑀Bi²𝐻²(1+Bi𝐻)³𝑑²𝐻𝑑𝑦²+𝑀Bi(1+Bi𝐻)³𝑑𝐻𝑑𝑦²−3𝑀Bi²𝐻(1+Bi𝐻)⁴𝑑𝐻𝑑𝑦²𝜂=𝜆−1−3𝑖𝜀𝑘ReΓ2𝐻²+5𝑖𝜀𝑘Re8𝐻⁴𝜂, (15) provided that |𝑒^𝜆𝑡|<(1/Bi|𝜁𝜂|)(1+Bi|𝐻|).

For the flow on a planar wall, the terms contributing to the real and imaginary parts of the complex eigenvalue 𝜆 are found as (𝐻=1, 𝑙=0, and 𝜂 is a constant)𝐶_𝑅=𝜀𝑘²65+Re+3ReΓ−cot𝛽𝑀Bi2(1+Bi)²−Re𝑊3𝑘²×1+𝜀²𝑘²3ReΓ2−5Re8⁻¹,𝐶_𝐿=3+𝜀3ReΓ2−5Re8𝐶_𝑅. (16a)(16b)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 (𝐶_𝑅>0) from (16a) becomes ()Re>cot𝛽−𝑀Bi/21+Bi²6/5+3Γ−𝑘²(),𝑊/3 (17) 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 𝑊 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 𝑦 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 𝑝(𝑦)(𝑑²𝜂/𝑑𝑦²)+𝑞(𝑦)𝜂=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 𝑦≠(1/2𝜋)cos⁻¹(1/𝐴_𝑙).

Since the coefficients of (15) are periodic with period 1, according to Floquet theory, 𝜂 can be expressed in the form [71] 𝜂(𝑦)=𝜗(𝑦)𝑒^𝑖𝜍𝑦, where 𝜗(𝑦) 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 𝜗(𝑦), when the topography is uneven, the eigenvalue problem (15) can be solved numerically to determine 𝐶_𝑅(𝑘,𝜍) and 𝐶_𝐿(𝑘,𝜍). Also, note that (15) remains invariant under the transformation 𝑦→−𝑦, suggesting that (15) is symmetric about the (𝑥,𝑧) plane (remark that (11) is actually symmetric about this plane), and hence, one can concentrate only in the region 𝑦∈[0,1/2] and demand 𝑑𝜂/𝑑𝑦=0, 𝑑³𝜂/𝑑𝑦³=0 at the end points. Another alternative would be the application of symmetry boundary conditions along 𝑦=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 𝑙=𝐴_𝑙cos(2𝜋𝑦) is defined and periodic in 𝑦 is mapped onto [−1,1] using the linear transformation 𝑌=2𝑦−1. The collocation points 𝑌_𝑗 in the interval [−1,1] are distributed according to the cosine law (where 𝑁 is a positive even integer): 𝑌_𝑗=cos𝑗𝜋𝑁,𝑗=0,1,2,…,𝑁−1,𝑁. (18) The derivatives 𝑑/𝑑𝑌 are replaced by a matrix 𝐷 of dimension (0∶𝑁,0∶𝑁), whose elements are𝐷()=𝑐𝑟,𝑠^𝑟(−1)^𝑟+𝑠𝑐_𝑠𝑌_𝑟−𝑌_𝑠𝐷()𝑌,0≤𝑟,𝑠≤𝑁,𝑟≠𝑠,𝑟,𝑟=−^𝑟21−𝑌²_𝑟𝐷()()=,1≤𝑟≤𝑁−1,0,0=−𝐷𝑁,𝑁2𝑁²+16, (19a)(19b)(19c)where 𝑌_𝑟,𝑌_𝑠 are distributed according to the cosine law (18). Further, 𝑐₀=𝑐_𝑁=2, 𝑐_𝑠=1, and 1≤𝑠≤𝑁−1. The vector obtained by multiplying 𝜂 with the matrix 𝐷 is the discrete approximation of 𝑑𝜂/𝑑𝑌. Higher order derivatives are approximated and obtained by repeated multiplication by the matrix 𝐷. For the cosine law based on the odd values of 𝑁 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 of 𝐏=(𝛽, 𝐴_𝑙, Re, 𝑊, Γ, 𝑀, 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 of 𝐶_𝑅 correspond to 𝜍=0, so that the most unstable modes are those for which 𝜗(𝑦) is periodic. Therefore, the dependence of 𝐶_𝑅 on 𝑘 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 𝐻(𝑦) theoretically constructed using the imposition 𝐴_𝑙≪1 can be checked and (ii) solutions at 𝐴_𝑙∼𝒪(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 𝐻(𝑦). It is observed that the numerical solution of 𝐻(𝑦) 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 𝑍(𝑦), 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, 𝑁≥40 was considered in the spectral technique, and 𝑁≥200 was chosen to implement the finite difference scheme.

(a)

(b)

(c)

Figure 2 

Free surface profiles corresponding to the primary flow on an uneven substrate with 𝐴𝑙=0.5, Re=5, 𝑀=5, and Bi=0.1: (a) cot𝛽=2, 𝑊=1; (b) cot𝛽=2, 𝑊=10, and (c) cot𝛽=0, 𝑊=1.

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 (Figure 3). 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 𝐻_max|_𝐴𝑙>0>𝐻_max|_𝐴𝑙=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 𝐴_𝑙, the variation in We shows either destabilizing or stabilizing tendency on the flow.

(a)

(b)

(c)

Figure 3 

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)𝑊=1, (2)𝑊=5, and (3)𝑊=10. In a(i)–c(i), solid lines represent 𝐴𝑙=0, and broken lines are the results for 𝐴𝑙=0.4, and (ii) and (iii) are for 𝐴𝑙=0.7 and 𝐴𝑙=0.9, respectively.

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 Figure 4 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 𝐴_𝑙=0. For example, referring to Figure 4(c) (ii) and Figure 4(c) (iii), it is seen that 𝐴_𝑙=0.4 yields larger growth rate values than 𝐴_𝑙=0, and the range of unstable wavenumbers in Figure 4(c) (ii) corresponding to 𝐴_𝑙=0 (unstable 𝑘∈[0,0.933]) is smaller than that for 𝐴_𝑙=0.4 (𝑘∈[0,0.9445]). For Figure 4(c) (iii), unstable 𝑘 for 𝐴_𝑙=0 lies in the interval [0,0.664], and that for 𝐴_𝑙=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.

(a)

(b)

(c)

Figure 4 

Dispersion of linear waves for isothermal Newtonian and viscoelastic films, temporal growth rate at Re=5, (a) Γ=0, (b) Γ=0.02 and (c) Γ=0.06; (1)𝐴𝑙=0, (2)𝐴𝑙=0.4, (3)𝐴𝑙=0.7, and (4)𝐴𝑙=0.9; (i) 𝑊=1, (ii) 𝑊=5, and (iii) 𝑊=10.

(a)

(b)

(c)

Figure 5 

Variation of growth rate with wavenumber for a vertically falling viscoelastic isothermal film with Γ=0.06 and Re=0.5, (1)𝑊=1, (2)𝑊=5, and (3)𝑊=10. In (i), solid lines represent 𝐴𝑙=0, and broken lines are the results for 𝐴𝑙=0.4, and (ii) and (iii) are for 𝐴𝑙=0.7 and 𝐴𝑙=0.9, respectively.

To substantiate the information availed from Figures 3, 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 𝐴_𝑙 when the viscoelastic effects increase at Re=5, as already discussed in relation to Figure 4.

(a)

(b)

Figure 6 

Stability chart in 𝑘−𝐴𝑙 plane for a vertically falling isothermal film: (a) Γ=0, (b) Γ=0.06; (i) Re=0.5, (ii) Re=1, and (1)𝑊=1, (2)𝑊=5, and (3)𝑊=10 10; S-Stable region, US-Unstable region.

(a)

(b)

Figure 7 

Neutral stability curves of an isothermal film at Re=5 and 𝛽=90∘: (a) 𝑊=5, (b) 𝑊=10, and (1)Γ=0, (2)Γ=0.02, and (3)Γ=0.06.

In Figure 8, the phase velocity curves varying with 𝑘 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 𝐶_𝑅>0, 𝐶_𝐿 would increase as a function of Γ only when Γ>5/12. Otherwise, despite 𝐶_𝑅>0, 𝐶_𝐿 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)𝐏(𝜋/2,0,5,5,0,0,0), (2)𝐏(𝜋/2,0.4,5,5,0,0,0), (3)𝐏(𝜋/2,0,5,5,0.02,0,0), (4)𝐏(𝜋/2,0.4,5,5,0.02,0,0), (5)𝐏(𝜋/2,0,5,5,0.06,0,0), and (6)𝐏(𝜋/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 𝐴_𝑙 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 𝐴_𝑙=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 𝐴_𝑙. However, it is also inferred from the figure that as 𝐴_𝑙 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 𝐴_𝑙=0 and 0.4 in Figure 11). The effect of thermocapillarity in this case is to increase the growth rate. However, at moderate 𝐴_𝑙=0.9 increase in heating stabilizes the flow, similar to the case discussed in Figure 9(d).

(a)

(b)

(c)

(d)

Figure 9 

The growth rate of the amplitudes at Re=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 𝑀=0, and the even numbers are for 𝑀=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 𝑀=0, it corresponds to an isothermal problem, and the effect of Bi is ignored.

Figure 10 

Influence of Marangoni stress on the neutral stability curves at Re=0.5 and Bi=0.1 on a vertically held substrate: (1)𝑊=5 and (2)𝑊=10. Broken lines are for 𝑀=0, and the solid lines represent 𝑀=10.

Figure 11 

Growth rate as a function of wavenumber at Re=5 and Bi=0.1.

Figure 12 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 𝐴_𝑙=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 𝐶_𝑅max=0.27 and 𝐶_𝑅max=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 𝑊 and 𝑀 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).

(a)

(b)

(c)

(d)

Figure 12 

Linear amplification rate: (1)𝐏(𝜋/2,𝐴𝑙,5,5,0,0,0), (2)𝐏(𝜋/2,𝐴𝑙,5,5,0.06,0,0), (3)𝐏(𝜋/2,𝐴𝑙,5,5,0,5,0.1), (4)𝐏(𝜋/2,𝐴𝑙,5,5,0.06,5,0.1), (5)𝐏(𝜋/2,𝐴𝑙,5,5,0,5,0.5) and (6)𝐏(𝜋/2,𝐴𝑙,5,5,0.06,5,0.5); (a) 𝐴𝑙=0, (b) 𝐴𝑙=0.4, and (c) 𝐴𝑙=0.9. In (d), the curves correspond to (2) and (6) from (a) and (b).

Figure 13 

Curves corresponding to the growth rate on a vertical topography with 𝐴𝑙=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 Figure 2 to get an insight into the physical mechanism involved. In this figure, it is observed that the free surface flattens up at 𝑧=1 for a large value of 𝑊 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 𝑢_0𝑧𝑧+3=0,𝑢₀,𝑢=0on𝑧=𝑙𝑥,𝑦_0𝑧=0on𝑧=1 (20) lead to 𝑢₀()+3𝑦,𝑧=3𝑧−𝑙2𝑙²−𝑧². (21) An integration of (21) yields the base flow rate as 𝑞=¹_𝑙𝑢₀𝑦,𝑧𝑑𝑧=(1−𝑙)³ (22) and the averaged flow rate as 𝑞_𝑎=¹₀𝑞𝑑𝑦=1+3𝐴²_𝑙2, (23) 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 𝑥-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 over 𝐴_𝑙, 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.

Appendices

A. Stress Tensor Components

The expressions for the components of the stress tensor 𝕋 in the explicit form are 𝜏_𝑥𝑥=−𝑝+2𝜇𝑢_𝑥−2Γ₀×𝑢_𝑥𝑡+𝑢𝑢_𝑥𝑥+𝑣𝑢_𝑥𝑦+𝑤𝑢_𝑥𝑧−2𝑢²_𝑥−𝑢²_𝑦−𝑢²_𝑧−𝑢_𝑦𝑣_𝑥−𝑢_𝑧𝑤_𝑥,𝜏_𝑦𝑦=−𝑝+2𝜇𝑣_𝑦−2Γ₀×𝑣_𝑦𝑡+𝑢𝑣_𝑥𝑦+𝑣𝑣_𝑦𝑦+𝑤𝑣_𝑦𝑧−𝑣²_𝑥−2𝑣²_𝑦−𝑣²_𝑧−𝑢_𝑦𝑣_𝑥−𝑣_𝑧𝑤_𝑦,𝜏_𝑧𝑧=−𝑝+2𝜇𝑤_𝑧−2Γ₀×𝑤_𝑧𝑡+𝑢𝑤_𝑥𝑧+𝑣𝑤_𝑦𝑧+𝑤𝑤_𝑧𝑧−𝑤²_𝑥−𝑤²_𝑦−2𝑤²_𝑧−𝑢_𝑧𝑤_𝑥−𝑣_𝑧𝑤_𝑦,𝜏_𝑥𝑦=𝜏_𝑦𝑥=𝜇𝑢_𝑦+𝑣_𝑥−Γ₀×𝑢_𝑦𝑡+𝑣_𝑥𝑡+𝑢𝑢_𝑥𝑦+𝑣_𝑥𝑥+𝑣𝑢_𝑦𝑦+𝑣_𝑥𝑦+𝑤𝑢_𝑦𝑧+𝑣_𝑥𝑧−3𝑢_𝑥𝑣_𝑥−3𝑢_𝑦𝑣_𝑦−2𝑢_𝑧𝑣_𝑧−𝑢_𝑥𝑢_𝑦−𝑣_𝑥𝑣_𝑦−𝑣_𝑧𝑤_𝑥−𝑢_𝑧𝑤_𝑦,𝜏_𝑥𝑧=𝜏_𝑧𝑥=𝜇𝑢_𝑧+𝑤_𝑥−Γ₀×𝑢_𝑧𝑡+𝑤_𝑥𝑡+𝑢𝑢_𝑥𝑧+𝑤_𝑥𝑥+𝑣𝑢_𝑦𝑧+𝑤_𝑥𝑦+𝑤𝑢_𝑧𝑧+𝑤_𝑥𝑧−3𝑢_𝑥𝑤_𝑥−3𝑢_𝑧𝑤_𝑧−2𝑢_𝑦𝑤_𝑦−𝑢_𝑥𝑢_𝑧−𝑤_𝑥𝑤_𝑧−𝑣_𝑥𝑤_𝑦−𝑢_𝑦𝑣_𝑧,𝜏_𝑦𝑧=𝜏_𝑧𝑦=𝜇𝑣_𝑧+𝑤_𝑦−Γ₀×𝑣_𝑧𝑡+𝑤_𝑦𝑡+𝑢𝑣_𝑥𝑧+𝑤_𝑥𝑦+𝑣𝑣_𝑦𝑧+𝑤_𝑦𝑦+𝑤𝑣_𝑧𝑧+𝑤_𝑦𝑧−3𝑣_𝑧𝑤_𝑧−2𝑣_𝑦𝑤_𝑦−2𝑣_𝑥𝑤_𝑥−𝑢_𝑦𝑤_𝑥−𝑤_𝑦𝑤_𝑧−𝑢_𝑧𝑣_𝑥−𝑣_𝑦𝑣_𝑧−𝑣_𝑦𝑤_𝑥. (A.1)

B. Dimensionless Equations and Solutions at the Leading and First Order

The nondimensional equations accurate up to 𝒪(𝜀) are given by 𝑢_𝑥+𝑣_𝑦+𝑤_𝑧𝑢=0,𝜀Re_𝑡+𝑢𝑢_𝑥+𝑣𝑢_𝑦+𝑤𝑢_𝑧=3−𝜀Re𝑝_𝑥+𝑢_𝑧𝑧×+𝜀ReΓ2𝑢_𝑧𝑢_𝑥𝑧+2𝑣_𝑧𝑢_𝑦𝑧+2𝑤_𝑧𝑢_𝑧𝑧−𝑢_𝑧𝑧𝑡+𝑢_𝑥𝑢_𝑧𝑧+𝑢_𝑦𝑣_𝑧𝑧+𝑢_𝑧𝑤_𝑧𝑧−𝑢𝑢_𝑥𝑧𝑧−𝑣𝑢_𝑦𝑧𝑧−𝑤𝑢_𝑧𝑧𝑧,𝜀Re𝑣_𝑡+𝑢𝑣_𝑥+𝑣𝑣_𝑦+𝑤𝑣_𝑧=−𝜀Re𝑝_𝑦+𝑣_𝑧𝑧×+𝜀ReΓ2𝑢_𝑧𝑣_𝑥𝑧+2𝑣_𝑧𝑣_𝑦𝑧+2𝑤_𝑧𝑣_𝑧𝑧−𝑣_𝑧𝑧𝑡+𝑣_𝑥𝑢_𝑧𝑧+𝑣_𝑦𝑣_𝑧𝑧+𝑣_𝑧𝑤_𝑧𝑧−𝑢𝑣_𝑥𝑧𝑧−𝑣𝑣_𝑦𝑧𝑧−𝑤𝑣_𝑧𝑧𝑧,𝜀𝑤_𝑧𝑧=3cot𝛽+Re𝑝_𝑧,𝑇_𝑧𝑧=𝜀Pe𝑇_𝑡+𝑢𝑇_𝑥+𝑣𝑇_𝑦+𝑤𝑇_𝑧. (B.1) The corresponding boundary conditions on 𝑧=𝑙(𝑥,𝑦) are 𝑢=0,𝑣=0,𝑤=0,𝑇=1, (B.2) and those on 𝑧=𝑓(𝑥,𝑦,𝑡) are ℎ_𝑡+𝑢𝑓_𝑥+𝑣𝑓_𝑦𝑇=𝑤,_𝑧+Bi𝑇=0,𝑝+2𝜀𝑓Re^𝑥𝑢_𝑧+𝑓_𝑦𝑣_𝑧−𝑤_𝑧=−𝜀²𝑊𝑒(1−Ca𝑇)𝑓_𝑥𝑥+𝑓_𝑦𝑦𝑢=0,_𝑧+𝜀𝑀𝑇_𝑥+𝑓_𝑥𝑇_𝑧×=𝜀ReΓ2𝑓_𝑥𝑢²_𝑧+2𝑓_𝑦𝑢_𝑧𝑣_𝑧+𝑢_𝑧𝑡+𝑢𝑢_𝑥𝑧+𝑣𝑢_𝑦𝑧+𝑤𝑢_𝑧𝑧−3𝑢_𝑧𝑤_𝑧−𝑢_𝑥𝑢_𝑧−𝑢_𝑦𝑣_𝑧,𝑣_𝑧+𝜀𝑀𝑇_𝑦+𝑓_𝑦𝑇_𝑧×=𝜀ReΓ2𝑓_𝑦𝑣²_𝑧+2𝑓_𝑥𝑢_𝑧𝑣_𝑧+𝑣_𝑧𝑡+𝑢𝑣_𝑥𝑧+𝑣𝑣_𝑦𝑧+𝑤𝑣_𝑧𝑧−3𝑣_𝑧𝑤_𝑧−𝑣_𝑥𝑢_𝑧−𝑣_𝑦𝑣_𝑧. (B.3)

Equations and BC’s at 𝒪(1) are𝑢_0𝑧𝑧+3=0,𝑣_0𝑧𝑧=0,Re𝑝_0𝑧𝑇+3cot𝛽=0,_0𝑧𝑧=0,𝑢_0𝑥+𝑣_0𝑦+𝑤_0𝑧𝑢=0,₀=𝑣₀=𝑤₀=0,𝑇₀,𝑢=1on𝑧=𝑙𝑥,𝑦_0𝑧=0,𝑣_0𝑧=0,𝑝₀=−𝜀²We𝑓_𝑥𝑥+𝑓_𝑦𝑦,𝑇on𝑧=𝑓𝑥,𝑦,𝑡_0𝑧+Bi𝑇₀=0on𝑧=𝑓𝑥,𝑦,𝑡. (B.4a)(B.4b)(B.4c)(B.4d)(B.4e)The solutions are𝑢₀=32𝑙²−𝑧²(),𝑣+3𝑓𝑧−𝑙₀𝑤=0,₀=3𝑓_𝑥2(𝑧−𝑙)²+3(𝑧−𝑙)ℎ𝑙_𝑥,𝑝₀=3cot𝛽−𝜀Re𝑓−𝑧²We𝑓_𝑥𝑥+𝑓_𝑦𝑦,𝑇₀=1+Bi𝑓−𝑧.1+Biℎ (B.5a)(B.5b)(B.5c)(B.5d)(B.5e)

Equations and BC’s at 𝒪(𝜀) areRe𝑢_0𝑡+𝑢₀𝑢_0𝑥+𝑣₀𝑢_0𝑦+𝑤₀𝑢_0𝑧=−Re𝑝_0𝑥+𝑢_1𝑧𝑧×+ReΓ2𝑢_0𝑧𝑢_0𝑥𝑧+2𝑣_0𝑧𝑢_0𝑦𝑧+2𝑤_0𝑧𝑢_0𝑧𝑧−𝑢_{0𝑧𝑧𝑡}+𝑢_0𝑥𝑢_0𝑧𝑧+𝑢_0𝑦𝑣_0𝑧𝑧+𝑢_0𝑧𝑤_0𝑧𝑧−𝑢₀𝑢_{0𝑥𝑧𝑧}−𝑣₀𝑢_{0𝑦𝑧𝑧}−𝑤₀𝑢_{0𝑧𝑧𝑧},Re𝑣_0𝑡+𝑢₀𝑣_0𝑥+𝑣₀𝑣_0𝑦+𝑤₀𝑣_0𝑧=−Re𝑝_0𝑦+𝑣_1𝑧𝑧×+ReΓ2𝑢_0𝑧𝑣_0𝑥𝑧+2𝑣_0𝑧𝑣_0𝑦𝑧+2𝑤_0𝑧𝑣_0𝑧𝑧−𝑣_{0𝑧𝑧𝑡}+𝑣_0𝑥𝑢_0𝑧𝑧+𝑣_0𝑦𝑣_0𝑧𝑧+𝑣_0𝑧𝑤_0𝑧𝑧−𝑢₀𝑣_{0𝑥𝑧𝑧}−𝑣₀𝑣_{0𝑦𝑧𝑧}−𝑤₀𝑣_{0𝑧𝑧𝑧},Re𝑝_1𝑧−𝑤_0𝑧𝑧𝑇=0,Pe_0𝑡+𝑢₀𝑇_0𝑥+𝑣₀𝑇_0𝑦+𝑤₀𝑇_0𝑧=𝑇_1𝑧𝑧𝑢_1𝑥,+𝑣_1𝑦+𝑤_1𝑧𝑢=0,₁=𝑣₁=𝑤₁=0,𝑇₁=0on𝑧=𝑙𝑥,𝑦,Re𝑝₁=2𝑤_0𝑧−𝑢_0𝑧𝑓_𝑥−𝑣_0𝑧𝑓_𝑦,𝑢on𝑧=𝑓𝑥,𝑦,𝑡_1𝑧+𝑀𝑇_0𝑥+𝑓_𝑥𝑇_0𝑧×=ReΓ2𝑢²_0𝑧𝑓_𝑥+2𝑢_0𝑧𝑣_0𝑧𝑓_𝑦+𝑢_0𝑧𝑡+𝑢₀𝑢_0𝑥𝑧+𝑣₀𝑢_0𝑦𝑧+𝑤₀𝑢_0𝑧𝑧−3𝑢_0𝑧𝑤_0𝑧−𝑢_0𝑥𝑢_0𝑧−𝑢_0𝑦𝑣_0𝑧,𝑣on𝑧=𝑓𝑥,𝑦,𝑡_1𝑧+𝑀𝑇_0𝑦+𝑓_𝑦𝑇_0𝑧×=ReΓ2𝑣²_0𝑧𝑓_𝑦+2𝑢_0𝑧𝑣_0𝑧𝑓_𝑥+𝑣_0𝑧𝑡+𝑢₀𝑣_0𝑥𝑧+𝑣₀𝑣_0𝑦𝑧+𝑤₀𝑣_0𝑧𝑧−3𝑣_0𝑧𝑤_0𝑧−𝑣_0𝑥𝑢_0𝑧−𝑣_0𝑦𝑣_0𝑧,𝑇on𝑧=𝑓𝑥,𝑦,𝑡_1𝑦+Bi𝑇₁=0on𝑧=𝑓𝑥,𝑦,𝑡. (B.6a)(B.6b)(B.6c)(B.6d)(B.6e)(B.6f)(B.6g)(B.6h)(B.6i)(B.6j)The solutions are𝑢₁=Re8×𝑓−𝑧−ℎ−4𝑓−𝑧²+8ℎ𝑓−𝑧+ℎℎ_𝑡+9𝑓𝑧ℎ2ℎ+𝑧−𝑓+3𝑓ℎ𝑓²+3ℎ²𝑧−3𝑓ℎ−3ℎ𝑧²+3𝑧ℎ+3ℎ²−3𝑓ℎ+9ℎ⁴ℎ_𝑥+32ReΓ𝑧−𝑓+ℎ3ℎ_𝑥𝑓ℎ−𝑧ℎ+3ℎ²+2ℎ_𝑡+𝑀Biℎ_𝑥(1+Biℎ)²+𝑧−𝑓+ℎ3cot𝛽𝑓_𝑥2−𝜀²WeRe2𝑓_𝑥𝑦𝑦+𝑓_𝑥𝑥𝑥×𝑧−𝑓²−ℎ²,𝑣₁1=−2𝜀²WeRe𝑓_𝑥𝑥𝑦+𝑓_𝑦𝑦𝑦−3cot𝛽𝑓_𝑦×𝑧−𝑓²−ℎ²+𝑀Biℎ_𝑦(1+Biℎ)²,𝑤𝑧−𝑓+ℎ₁=𝑢_1𝑥+𝑣_1𝑦𝑝𝑑𝑧,₁3=−𝑓Re𝑧−𝑓_𝑥−6Reℎℎ_𝑥,𝑇₁1=−120BiPeℎ_𝑡(1+Biℎ)³×−60𝑧²ℎ²+40𝑓ℎ³+20𝑓³ℎ−40𝑧ℎ³−60𝑓²ℎ²−20𝑧³ℎ−60𝑧𝑓²ℎ+120𝑧𝑓ℎ²+60𝑧²𝑓ℎBi²+−60𝑧²ℎ+60𝑧²𝑓−60𝑓²ℎ−20𝑧³+40ℎ³+20𝑓³+120𝑧𝑓ℎ−60𝑧𝑓²−1Bi120BiPeℎ_𝑥(1+Biℎ)³×−90𝑧³ℎ³−270𝑧𝑓²ℎ³+45𝑧𝑓⁴ℎ+270𝑧²ℎ³𝑓+90𝑓³ℎ³−180𝑧²ℎ⁴−99𝑧ℎ⁵−90𝑧²𝑓³ℎ−180𝑓²ℎ⁴+9𝑧⁵ℎ−9𝑓⁵ℎ+360𝑧𝑓ℎ⁴+99𝑓ℎ⁵−45𝑧⁴𝑓ℎ+90𝑧³𝑓²ℎBi²+60𝑓𝑧³ℎ−270𝑧²ℎ³−15𝑓⁴ℎ−45𝑧ℎ⁴−150𝑧³ℎ²+450𝑧²ℎ²𝑓+9𝑧⁵−9𝑓⁵+99ℎ⁵+45𝑓ℎ⁴−450𝑧𝑓²ℎ²+90𝑧³𝑓²+150𝑓³ℎ²−90𝑧²𝑓²ℎ−270𝑓²ℎ³−90𝑧²𝑓³+540𝑧𝑓ℎ³−15𝑧⁴ℎ−45𝑧⁴𝑓+45𝑧𝑓⁴+60𝑧𝑓³ℎBi+180𝑧𝑓ℎ²+60𝑓³ℎ−180𝑧𝑓²ℎ+180𝑧²𝑓ℎ−90𝑓²ℎ²+60𝑓𝑧³+45ℎ⁴+60𝑧𝑓³−90𝑧²ℎ²−90𝑧²𝑓²−15𝑧⁴−15𝑓⁴−60𝑧³ℎ. (B.7a)(B.7b)(B.7c)(B.7d)(B.7e)