Computational-Analysis of the Non-Isothermal Dynamics of the Gravity-Driven Flow of Viscoelastic-Fluid-Based Nanofluids Down an Inclined Plane

Nomenclature

1  Introduction

Dispersion of thermally conductive, tiny metallic particles within a fluid (such as water, oil, etc.) is the most obvious and effective method of improving/enhancing the heat-transfer-rate (HTR) characteristics and thermal-conductivity properties of the fluid. The size and texture of the metallic particles are clearly fundamental; large particles can lead to sedimentation, clogging, etc., while coarse-grained particles can cause abrasion, etc. For these reasons, nanometer-size particles (nanoparticles) are used [1–24].

As in the cited references, the term nanofluid is herein used to describe the suspension of solid nanoparticles in a base fluid. A variety of industrial applications (say heating and cooling) and medical applications (chemotherapy, etc.) of nanofluids are well documented in the literature and summarized by the cited references. Noting the fundamental importance of the rheology of the base fluid to the modern applications of nanofluids, the studies, say in [1–3] focus attention on non-Newtonian (specifically polymeric) base-fluids under various flow conditions. The widespread importance of non-Newtonian fluids, both the Generalized Newtonian Fluids (GNF) as well as the polymeric (viscoelastic) fluids, in contemporary application is quite straightforward for example [1–3,6,8,25–34].

The current study builds on the investigations in [1–3] and extends these investigations to free-surface, gravity-driven, thin-film flows with convective heat exchange at the free-surface. Additionally and alternatively, the current work extends the combined (particle-free) fluid-dynamics studies in [25] and [29] to the inclusion of solid nanoparticles and hence raising the cited works to nanofluid-dynamics.

The current work employs solution methodologies based on semi-implicit numerical-schemes derived from the finite difference methods (FDM). The computational solution algorithms are implemented in the MATLAB software. It is noted that other alternative numerical methodologies have been utilized for the computational solution processes of similar nanofluid flow problems. Such alternative numerical methodologies have been applied for hybrid nanofluid flow [35,36], for fractional viscoelastic fluid [37,38], and for fractional viscoelastic hybrid nanofluid [39]. For example, the work in [40] investigates a Generalized Newtonian fluid (the Casson fluid) based nanofluid subjected to thermal radiation, Joule heating, and magnetic effects flowing over an inclined porous stretching sheet. Similar work in [41] investigates thermal and velocity slip effects on Casson nanofluid flow over an inclined permeable stretching cylinder via collocation methods. The work in [42] gives an extensive review of strategies based, among others, on nanofluids for solving “internal heat generation problem” for the optimization of heat transfer in electronic devices. The work in [43] explores the significant potential of nanofluids, and specifically the influence of nano-particle size and distribution, in solar energy harvesting. The importance of viscoelastic fluid in heat generation/absorption is also explored in [44] using Maxwell fluids.

The following sequence is adopted in the paper. Section 2 outlines the description of the physical and mathematical models. The development and implementation of the numerical and computational algorithms as well as the real efficacy, accuracy, and convergence of the computational methodologies are given in Section 3. The main results are presented graphically and discussed qualitatively in Section 4. Concluding remarks follow in Section 5.

2  Problem Formulation

A schematic of the model problem is sketched in Fig. 1. A thin film of nanofluid (VFBN) of height h∗ (where the height is in the direction normal to the inclined wall) flows down an inclined plane. The inclined wall/plane makes an angle θ with the horizontal.

Figure 1: Schematic of the model problem

The x∗ -axis is taken parallel to the inclined wall and the y∗ -axis is taken perpendicular to this incline. The superscript ( ∗ ) denotes dimensional variables. The variables and parameters are otherwise similar to those described in [1–3].

The solid boundary (the inclined wall) is kept at a constant temperature Tw∗ and convective cooling is assumed at the free-surface, governed by Newton’s law of cooling.

The motion of the fluid is exclusively gravity-driven and hence the pressure gradient in the x∗ -direction stays zero. Alternatively the pressure in the x∗ -direction stays constant. Naturally, the pressure in the y∗ -direction is non-constant increasing from atmospheric pressure at the free surface to the maximum fluid pressure at the bottom, i.e., at the inclined wall.

No-slip velocity boundary conditions are assumed along the rigid inclined wall. The velocity boundary conditions at the free-surface naturally also arise from the zero-shear-rate requirements.

Following the notations of [1–3,29], the governing PDEs (in dimensional form) for the VFBN are,

∇∗⋅u∗=0, (1)

ρnf∗Du∗Dt∗=−∇∗p∗+∇∗⋅(σ∗)+ρnf∗g∗F, (2)

(ρcp)nf∗DT∗Dt∗=−∇∗⋅ϕq∗+QD∗+r∗. (3)

Here F=(sin⁡θ,−cos⁡θ) denotes the body force due to gravity and g∗ is the acceleration due to gravity. The heat-source terms, denoted by r∗ , will subsequently be neglected in this study as their effects have been comprehensively investigated in previous studies, as in [1].

Dimensionless parameters

The equations are studied in non-dimensional form via the following dimensionless-parameters; Deborah-number (De), Reynolds-number (Re), Prandtl-number (Pr), Brinkman-number (Br), Peclet-number (Pe = Re ⋅ Pr), activation-energy parameter ( α ), Grashof-number (Gr), Biot-number (Bi), and ( β ) which represents the proportion of the polymer-viscosity compared to the fluid’s total-viscosity. These are defined as follows:

β=ηp∞∗η∞∗,De=λ∞∗U∞∗h∗,Re=ρf∗h∗U∞∗η∞∗,Pr=η∞cpf∗∗κf∗,Br=U∞∗2η∞∗κf∗αTW∗,Renf=ρnfρfRe,Penf=(ρcp)nf(ρcp)fPe,α=R∗TW∗E∗,Gr=g∗h∗U∞∗,Bi=H∗h∗k0∗. (4)

The subscript ( nf ) depicts nanofluid and the appropriate (nanofluid) quantity is obtained linearly from contributions of the base-fluid volume-fractions ( f ) and the nanoparticles ( s ) volume-fractions,

ρnf=φρs+(1−φ)ρf,(cpρ)nf=φ(cpρ)s+(1−φ)(cp,ρ)f,

where, φ is the volume fraction function. The resultant dimensionless governing equations are,

∇⋅u=0, (5)

RenfDuDt=−Renf∇p+∇⋅(σ)+RenfGrF, (6)

PenfDTDt=−∇⋅(κnf∇T)+BrQD, (7)

τ+ετ2+λ¯(T)De[ τ∇ − τDDt(ln(αT+1)) ]=ηp(T)(1−φ)5/2S.(8)

The expression for the appropriate mechanical dissipation term is,

QD=γτ:S+(1−γ)ηs(T)(1−φ)5S:S. (9)

As given in [1], the expressions for viscosities, relaxation-times, and thermal-conductivities, under non-isothermal conditions are,

ηs(T)=(1−β)e(−αT),(ηs)nf=ηs(T)(1−φ)5/2, (10)

ηp(T)=βe(−αT),(ηp)nf=ηp(T)(1−φ)5/2, (11)

η=ηp(T)+ηs(T),ηnf=η(1−φ)5/2, (12)

λ¯(T)=1αT+1e(−αT), (13)

κnf=κs−(ℵ−1)φ(κf−κs)+(1−ℵ)κfκs+φ(κf−κs)+(1−ℵ)κf(1+αA2T). (14)

2.1 Boundary and Initial Conditions

The start-up conditions (at time t = 0) and the wall-conditions (at the inclined wall, y = 0, and at the free-surface, y = 1) are,

T(0,y)=0,u(0,y)=0,τ(0,y)=0,0≤y≤1, (15)

u(t,0)=0,T(t,0)=0,∂∂yu(t,1)=1,∂∂yT(t,1)=−BiT(t,1),t≥0. (16)

Due to the hyperbolic nature of the equations for the polymeric stresses, the relevant boundary-conditions for these equations are reconstructed from the main flow [1–3,29].

3  Numerical Solution

The numerical-scheme adopted for the velocity component is derived from a semi-implicit FDM approach,

Renfu(n+1)−u(n)Δt=(ηs)nf(n)∂2∂y2u(n+ξ)+∂∂yτ12(n)+∂∂yu(n)∂∂y(ηs)nf(n)+RenfGrsin⁡θ, (17)

where

u(n+ξ)=(1−ξ)u(n)+ξu(n+1).

The velocity therefore updates at the new time-level, u(n+1) , via the recursive algorithm,

−r1uj−1(n+1)+(Renf+2r1)uj(n+1)−r1uj+1(n+1)=Renfuj(n)+(1−ξ)Δt(ηs)nf(n)∂2∂y2u(n)+Δt((1−β)∂∂y(ηs)nf(n)∂∂yu(n)+∂∂yτ12(n)+ΔtRenfGrsin⁡θ), (18)

where

r1=ξ(ηs)nf(n)ΔtΔy2.

This represents a (diagonally-dominant) tri-diagonal linear algebraic system of equations. The discretized temperature equation is obtained similarly,

Penf∂T∂t=κnf(n)∂2∂y2T(n+ξ)+∂∂yT(n)∂∂yκnf(n)+BrQD(n). (19)

The temperature therefore updates at the new time-level, T(n+1) , via the recursive algorithm,

−r2Tj−1(n+1)+(Penf+2r2)Tj(n+1)−r2Tj+1(n+1)=PenfTj(n)+(1−ξ)Δtκnf∂2∂y2T(n)+Δt∂∂yT(n)∂∂yκnf(n)+2(1−γ)ΔtBr(ηs)nf(n)(∂∂yu(n))2+2ΔtBrγτ12(n)∂∂yu(n), (20)

where

r2=ξκnf(n)ΔtΔy2.

The semi-implicit numerical scheme for the polymeric-stress, τ is,

τ(n+ξ)+ε(τ2)(n)+λ¯(n)Deτ(n+1)−τ(n)Δt=explicitterms.

The solutions for the tensor components, τ11(n+1) , τ12(n+1) , and τ22(n+1) therefore follow directly from algebraic manipulation,

(Deλ¯(n)+ξΔt)τ(n+1)=explicit‐terms. (21)

The explicit-terms for τ11 , τ12 and τ22 are respectively,

[λ¯(n)De−(1−ξ)Δt]τ11(n)+Δtλ¯(n)De[τ12(n)∂∂yu(n)+τ11(n)∂∂ylog⁡(1+αT(n))]−εΔt(τ112+τ122), (22)

[λ¯(n)De−(1−ξ)Δt]τ12(n)+Δtλ¯(n)De[τ22(n)∂∂yu(n)+τ12(n)∂∂ylog⁡(1+αT(n))]+Δt(ηp)nf∂∂yu(n)−εΔt(τ11τ12+τ12τ22), (23)

[λ¯(n)De−(1−ξ)Δt]τ22(n)+Δtλ¯(n)Deτ22(n)∂∂ylog⁡(1+αT(n))−εΔt(τ122+τ222). (24)

4  Results

Graphical results are presented for the VFBN temperature field (T), velocity field (u), and polymeric-stress components ( τ11,τ12andτ22 ), using the below list of default values,

α=0.01,Br=1,Prnf=1,Renf=1,Bi=1,Gr=1,De=2,β=0.2,γ=0.5,Δy=0.01,Δt=0.1,θ=45,t=50,ζ=1,φ=0.04,A2=0.2,ε=1,ℵ=3. (25)

4.1 Time Devolvement of Steady Solutions

The time development of flow variables from the initial states until steady solutions are reached, as is illustrated in Figs. 2 and 3. It can be noted that all solutions stably and ultimately settle to consistent steady states.

Figure 2: Time-development of profiles to steady-state with Δt=0.05

Figure 3: Time development of maximum flow quantities to steady-state

4.2 Time-Step and Mesh Size Convergence

Figs. 4–7 illustrate that the computational and numerical algorithms are independent of both time-step size and mesh size as required. The computational and numerical scheme specifically, and efficiently, reproduce the required steady-solutions for a wide range of expected time-step sizes and mesh sizes.

Figure 4: Time-step independence of maximum solutions

Figure 5: Time-step independence of solutions

Figure 6: Mesh-size independence of maximum solutions

Figure 7: Mesh-size independence of solutions

4.3 Code Validation

A similar investigation to the current one was conducted in [29] using the Oldroyd-B constitutive model and in the absence of nanoparticles. The Oldroyd-B model is obtained from the Giesekus model by taking the nonlinear parameters as zero, i.e., ε=0 . The absence of nanoparticles in [29] mathematically reduces to taking φ=0 . The study in [29] therefore used the ordinary Oldroyd-B constitutive model with constant thermal-conductivity, κnf≡κf , i.e., κs=ℵ=A2=0 . By taking κs=0 , ℵ=0 , ε=0 , φ=0 , and A2=0 in the current model, the VFBN results reduce to those for a normal viscoelastic (Oldroyd-B) fluid and are exactly the same as those in [29].

4.4 Sensitivity of Solutions to Embedded Parameters

A representative sample of the flow behavior of field variables with variations in the fluid and flow parameters is now presented. The sensitivity of the VFBN thermal conductivity to changes in φ is illustrated in Fig. 8. As expected, the VFBN thermal-conductivity increases with increasing nanofluid volume-fraction, φ .

Figure 8: Variation of VFBN thermal conductivity with φ

It therefore naturally follows that the VFBN temperature increases with increasing φ . This expected response of the VFBN temperature to variations in φ is illustrated in Fig. 9. An increase in fluid temperature directly results in a reduction in fluid viscosity and hence an increase in the fluid velocity. It is therefore expected that the VFBN velocity would increase with increasing VFBN temperature. This expected increase of the VFBN velocity to variations in φ is also illustrated in Fig. 9. Lastly, Fig. 9 demonstrates an increase in the polymeric-stresses with increasing nanofluid volume-fraction, φ .

Figure 9: Variation of VFBN flow quantities with φ

Similar responses (as with φ ) of the VFBN thermal-conductivity to variations in both the thermal-conductivity parameter, A2 , and the activation-energy parameter, α , are respectively illustrated in Figs. 10 and 11.

Figure 10: Variation of VFBN thermal conductivity with A2

Figure 11: Variation of VFBN thermal conductivity with α

Fig. 12 shows the response of flow quantities to variations in the nanofluid Prandtl number, Prnf . Noting that Prnf is inversely proportional to the strength of the heat sources, the VFBN thermal conductivity (subsequently therefore the VFBN temperature) should decrease as Prnf increases. These expected results illustrated in Figs. 12 and 13, respectively.

Figure 12: Effects of Prnf on VFBN thermal conductivity

Figure 13: Effects of Prandtl Number, Prnf , on VFBN temperature

The behaviour of VFBN thermal conductivity with variations in the Deborah number, De, is illustrated in Fig. 14. It is observed that the VFBN thermal-conductivity increases with increasing De.

Figure 14: Effects of De on VFBN thermal conductivity

4.5 Wall Shear-Stress and Wall Heat Transfer Rate

As was already observed in the previous section, the nanofluid (VFBN) velocity and temperature both increase with increasing φ,Prnf , and De, see Figs. 8, 9, 12–14. The wall shear-stress and wall heat transfer rate (at the inclined wall) follow the same patterns, they both also increase with increasing φ,Prnf , and De.

5  Concluding Remarks

Efficient and robust computational and numerical schemes based on semi-implicit FDM were employed to investigate the thermo-physical and fluid-dynamical characteristics of a non-isothermal, free-surface, thin-film, gravity-driven flow of a VFBN. The VFBN flows down an inclined-plane and is subjected to convective-cooling at the fluid-air interface, i.e., the free-surface. The numerical and computational algorithms were checked for convergence in both space and time and were also positively validated against the results in the existing literature. A single-phase nanofluid model, in which metallic nanoparticles of spherical shape are homogenously mixed (ensuring non-sedimentation) to a viscoelastic base fluid of the Giesekus type, was adopted. The results illustrate that the volume-fraction of the embedded nanoparticles play fundamental roles in the fluid-dynamical and thermodynamical properties of the VFBN. Specifically, the VFBN thermal-conductivity, VFBN temperature, VFBN velocity, and VFBN polymeric stresses (in particular, the first normal stress difference) all increase as the nanoparticle volume-fraction increases. The results are of fundamental industrial significance and importance, especially to heating and cooling applications where heat-transfer-rate (HTR) enhancement and thermal-conductivity improvement may be achieved via nanoparticles and nanofluidics.

Funding Statement: The authors received no specific funding for this study.

Conflicts of Interests: The authors declare that they have no conflicts of interest to report regarding the present study.