Computational Analysis of Shear Banding in Simple Shear Flow of Viscoelastic Fluid-Based Nanofluids Subject to Exothermic Reactions

Abstract

We investigated the shear banding phenomena in the non-isothermal simple-shear flow of a viscoelastic-fluid-based nanofluid (VFBN) subject to exothermic reactions. The polymeric (viscoelastic) behavior of the VFBN was modeled via the Giesekus constitutive equation, with appropriate adjustments to incorporate both the non-isothermal and nanoparticle effects. Nahme-type laws were employed to describe the temperature dependence of the VFBN viscosities and relaxation times. The Arrhenius theory was used for the modeling and incorporation of exothermic reactions. The VFBN was modeled as a single-phase homogeneous-mixture and, hence, the effects of the nanoparticles were based on the volume fraction parameter. Efficient numerical schemes based on semi-implicit finite-difference-methods were employed in MATLAB for the computational solution of the governing systems of partial differential equations. The fundamental fluid-dynamical and thermodynamical phenomena, such as shear banding, thermal runaway, and heat transfer rate (HTR) enhancement, were explored under relevant conditions. Important novel results of industrial significance were observed and demonstrated. Firstly, under shear banding conditions of the Giesekus-type VFBN model, we observed remarkable HTR and Therm-C enhancement in the VFBN as compared to, say, NFBN. Specifically, the results demonstrate that the VFBN are less susceptible to thermal runaway than are NFBN. Additionally, the results illustrate that the reduced susceptibility of the Giesekus-type VFBN to the thermal runaway phenomena is further enhanced under shear banding conditions, in particular when the nanofluid becomes increasingly polymeric. Increased polymer viscosity is used as the most direct proxy for measuring the increase in the polymeric nature of the fluid.

1. Introduction

Research into the enhancement of heat transfer rate (HTR) and thermal conductivity (Therm-C) characteristics of nanofluids has gained wide contemporary traction, see for example [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. The evidence of the efficacy of nanofluids (formed from mixtures of metallic nanoparticles and base fluids, such as water or oil) in HTR and Therm-C enhancement is now widespread and commonly accepted; hence, such nanofluids are now extensively used for these purposes.

Similar investigations have been conducted on the use of various (natural and synthetic) fluid properties for the enhancement of certain characteristics that are desirable for applications. Fluid viscoelasticity has indeed received wide attention in this regard, specifically for the improvement of certain industrial, domestic, and medical applications. Such widespread use of the viscoelastic (polymeric) fluids, for the improvement and enhancement of, say, HTR and Therm-C, emulsion polymerization, drag reduction, etc., can be found, for example, in [1,4,6,23,24,25,26,27,28,29,30,31,32,33,34].

The present study extends the investigations in [1] to the conditions of shear banding. Specifically, the current investigation explores the novel contributions of shear banding phenomena to HTR and Therm-C enhancement. Alternatively (and additionally), the current investigation extends the studies in [25,26,27] (all of which focused on the flow behavior and thermodynamics of particle-free viscoelastic fluids) to nanofluids, specifically VFBN modeled by the Giesekus constitutive model. The study in [25] investigates shear banding phenomena using the two-fluid Rolie-Poly viscoelastic constitutive model and the investigations in [26,27] employ the Johnson-Segalman viscoelastic constitutive model.

Shear banding phenomena in the shear flow of viscoelastic fluids represents observable and physical discontinuities in the shear rate profiles of the flow velocity and finds wide industrial and domestic application, say, in emulsion polymerization and drop break-up [34]. The main mechanism for shear banding in [25] falls under the category of what are termed flow-induced inhomogeneities. On the other hand, the mechanisms for shear banding in [26,27] fall under the category of constitutive instabilities. However, the researchers in [25] demonstrated that their finite volume method (FVM)-based numerical and computational algorithms (implemented on the OpenFOAM software) were capable of reproducing the shear banding characteristics via either mechanism, using the Johnson-Segalman, Giesekus, and Rolie-Poly viscoelastic constitutive models. Shear banding phenomena has been observed in the shear flow of polymeric fluids of various types, including, for example, worm-like micelle solutions [35,36], polymer solutions [37], foams [38], telechelic polymers [39], granular flows [40], soft glasses [41], and polymer melts [42].

The paper is structured as follows. Section 2 presents the physical and mathematical model. Section 3 outlines the numerical and computational algorithms. Section 4 lists the graphical and qualitative results. Section 5 furnishes the concluding remarks.

2. Problem Formulation

We consider plane Couette flow (simple shear flow) as depicted in the schematic in Figure 1.

The variables appearing in Figure 1 are the velocity field, $u(t,\mathit{x})$, the temperature field, $T(t,\mathit{x})$, the time, t, and the 2D Cartesian space coordinates, $\mathit{x}=(x,y)$.

2.1. Model Assumptions

• We assume the flow of a VFBN in a channel of infinite longitudinal extent. We therefore assume that the flow is fully developed in the x-direction and, hence, that flow quantities are independent of x.
  This allows us to focus our attention on the primary effects of shear banding on HTR and Therm-C enhancement without the complications of 2D (or indeed 3D) computations.
• We assume that the shear banding is driven by constitutive instabilities via the Giesekus viscoelastic constitutive model.
  The exact mechanisms of shear banding, whether via constitutive instabilities or via flow inhomogeneities, are still areas of active research. Indeed, even for shear banding via constitutive instabilities, at least two viscoelastic constitutive models have been advanced. None of these considerations however detract from the primary aim to investigate the broader effects of shear banding on HTR and Therm-C enhancement.
• We assume spherical nanoparticles that are homogeneously mixed with the base-fluid.
  The size, shape, distribution, orientation, etc., of the nanoparticles are still wide open areas with regard to investigating the optimal conditions for HTR and Therm-C enhancement. However, these considerations do not detract from the primary aim—to investigate the broader effects of nanoparticles on HTR and Therm-C enhancement.

2.2. Dimensionless Governing Equations

The development of the governing equations for the VFBN is adapted from [1]. Using the notation of [1], the relevant equations in dimensionless form are,

$$\begin{array}{cc}\hfill & \nabla ·\mathit{u}=0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & {\mathrm{Re}}_{nf}\frac{D\mathit{u}}{Dt}=-{\mathrm{Re}}_{nf}\nabla p+\nabla ·\mathit{\sigma },\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & {\mathrm{Pe}}_{nf}\frac{DT}{Dt}=-\nabla ·\left({\kappa }_{nf}\nabla T\right)+\mathrm{Br}{Q}_D+{\delta }_{1}exp\left(\frac{T}{1+\alpha T}\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & \mathit{\tau }+\epsilon {\mathit{\tau }}^2+\mathrm{De}\overline{\lambda }\left(T\right)\left[\stackrel{\nabla }{\mathit{\tau }}-\mathit{\tau }\frac{D}{Dt}\left(ln\left(1+\alpha T\right)\right)\right]=\frac{{\eta }_p\left(T\right)}{{\left(\sqrt{1-\varphi }\right)}^5}\mathit{S},\hfill \end{array}$$

(4)

with,

$$\begin{array}{cc}\hfill & Q_D=\gamma \mathit{S}:\mathit{\tau }+(1-\gamma )\frac{{\mathit{\eta }}_s\left(T\right)}{{\left(\sqrt{1-\varphi }\right)}^5}\mathit{S}:\mathit{S},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\eta }_s\left(T\right)=\left(1-\beta \right)exp(-\alpha T),{\left({\eta }_s\right)}_{nf}=\frac{{\eta }_s\left(T\right)}{{\left(\sqrt{1-\varphi }\right)}^5},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & {\eta }_p\left(T\right)=\beta e^{(-\alpha T)},{\left({\eta }_p\right)}_{nf}=\frac{{\eta }_p\left(T\right)}{{\left(\sqrt{1-\varphi }\right)}^5},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \eta ={\eta }_s\left(T\right)+{\eta }_p\left(T\right),{\eta }_{nf}=\frac{\eta }{{\left(\sqrt{1-\varphi }\right)}^5},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \overline{\lambda }\left(T\right)=\frac{1}{1+\alpha T}exp(-\alpha T),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\kappa }_{nf}=\frac{{\kappa }_s+(1-\chi ){\kappa }_f-(\chi -1)\varphi ({\kappa }_f-{\kappa }_s)}{{\kappa }_s+(1-\chi ){\kappa }_f+\varphi ({\kappa }_f-{\kappa }_s)}(1+\alpha A_{2}T),\hfill \end{array}$$

(10)

where,

$$\mathit{\sigma }=\mathit{\tau }+{\left({\eta }_s\right)}_{nf}\mathit{S},\mathit{S}=\frac{1}{2}\left(\nabla \mathit{u}+{(\nabla \mathit{u})}^T\right).$$

(11)

The non-dimensional variables and dimensionless parameters used in Equations (1)–(11) are defined in the nomenclature at the end of the article. In particular, p is the pressure field; $\epsilon$ is the Giesekus non-linear parameter; $\sigma$ is the total stress tensor; $S$ is the rate of deformation tensor; $\mathit{\tau }$ is the polymer stress tensor; ${\left({\eta }_s\right)}_{nf}$ is the solvent viscosity for the nanofluid; ${\left({\eta }_p\right)}_{nf}$ is the polymer viscosity for the nanofluid; ${\eta }_{nf}$ is the total viscosity for the nanofluid; ${\kappa }_{nf}$ is the thermal conductivity for the nanofluid; $\alpha$ is the activation energy parameter; $\beta$ is the polymer to total viscosity ratio; Br is the Brinkman number; ${\delta }_1$ is the Frank-Kamenetskii parameter; De is the Deborah number; Pr is the Prandtl number; Pe is the Peclet number; Re is the Reynolds number, and the subscript ${\left(\right)}_{nf}$ represent nanofluid. Each os the nanofluid quantities, ${\left(\right)}_{nf}$, is obtained as a linear combination of the contributions from the solid nanoparticles, ${\left(\right)}_s$, and the base-fluid, ${\left(\right)}_f$;

$${\left(\right)}_{nf}=\varphi {\left(\right)}_s+(1-\varphi ){\left(\right)}_f.$$

The notations ∇ and $D/Dt$ represent the usual gradient and material-derivative operators respectively.

2.3. Initial and Boundary Conditions

The following initial and boundary conditions apply

$$\begin{array}{cc}\hfill & u(0,y)=0,T(0,y)=0,\mathit{\tau }(0,y)=0,\mathrm{for}0\le y\le 1,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & u(t,0)=0,u(t,1)=1,T(t,0)=0,T(t,1)=0,\mathrm{for}t\ge 0.\hfill \end{array}$$

(13)

3. Numerical and Computational Algorithms

The semi-implicit numerical scheme for the longitudinal velocity component is

$$\begin{array}{c}\hfill {\mathrm{Re}}_{nf}\frac{u^{(n+1)}-u^{\left(n\right)}}{\Delta t}={\mathrm{Re}}_{nf}G+\frac{\partial }{\partial y}{\tau }_{12}^{\left(n\right)}+\left({\left({\eta }_s\right)}_{nf}^{\left(n\right)}\frac{{\partial }^2}{\partial y^2}{u}^{(n+\xi )}+\frac{\partial }{\partial y}{u}^{\left(n\right)}\frac{\partial }{\partial y}{\left({\eta }_s\right)}_{nf}^{\left(n\right)}\right),\end{array}$$

(14)

where

$$u^{(n+\xi )}=\xi u^{(n+1)}+(1-\xi )u^{\left(n\right)}.$$

Therefore, the longitudinal velocity updates in time to $u^{(n+1)}$, via

$$\begin{array}{cc}\hfill -r_1{u}_{j-1}^{(n+1)}+({\mathrm{Re}}_{nf}+2r_1)u_j^{(n+1)}-r_1{u}_{j+1}^{(n+1)}& ={\mathrm{Re}}_{nf}G+{\mathrm{Re}}_{nf}{u}_j^{\left(n\right)}+\Delta t\frac{\partial }{\partial y}{\tau }_{12}^{\left(n\right)}\hfill \\ \hfill & +(1-\xi )\Delta t{\left({\eta }_s\right)}_{nf}^{\left(n\right)}\frac{{\partial }^2}{\partial y^2}{u}^{\left(n\right)}\hfill \\ \hfill & +(1-\beta )\Delta t\frac{\partial }{\partial y}{u}^{\left(n\right)}\frac{\partial }{\partial y}{\left({\eta }_s\right)}_{nf}^{\left(n\right)},\hfill \end{array}$$

(15)

where

$$r_1=\xi {\left({\eta }_s\right)}_{nf}^{\left(n\right)}\frac{\Delta t}{\Delta y^2}.$$

A similar semi-implicit numerical scheme for the temperature is

$$\begin{array}{c}\hfill {\mathrm{Pe}}_{nf}\frac{\partial T}{\partial t}={\kappa }_{nf}^{\left(n\right)}\frac{{\partial }^2}{\partial y^2}{T}^{(n+\xi )}+\frac{\partial }{\partial y}{T}^{\left(n\right)}\frac{\partial }{\partial y}{\kappa }_{nf}^{\left(n\right)}+\mathrm{Br}{Q}_D^{\left(n\right)}+{\delta }_{1}exp\left(\frac{T^{\left(n\right)}}{1+\alpha T^{\left(n\right)}}\right).\end{array}$$

(16)

leading to,

$$\begin{array}{c}\hfill -r_2{T}_{j-1}^{(n+1)}+({\mathrm{Pe}}_{nf}+2r_2)T_j^{(n+1)}-r_2{T}_{j+1}^{(n+1)}={\mathrm{Pe}}_{nf}{T}_j^{\left(n\right)}+(1-\xi )\Delta t{\kappa }_{nf}\frac{{\partial }^2}{\partial y^2}{T}^{\left(n\right)}\\ \hfill +\Delta t\frac{\partial }{\partial y}{T}^{\left(n\right)}\frac{\partial }{\partial y}{\kappa }_{nf}^{\left(n\right)}+2(1-\gamma )\Delta t\mathrm{Br}{\left({\eta }_s\right)}_{nf}^{\left(n\right)}{\left(\frac{\partial }{\partial y}{u}^{\left(n\right)}\right)}^2+2\Delta t\mathrm{Br}\gamma {\tau }_{12}^{\left(n\right)}\frac{\partial }{\partial y}{u}^{\left(n\right)},\end{array}$$

(17)

where

$$r_2=\xi {\kappa }_{nf}\frac{\Delta t}{\Delta y^2}.$$

Finally, the semi-implicit numerical scheme for the polymeric-stress is

$${\mathit{\tau }}^{(n+\xi )}+\epsilon {\left({\mathit{\tau }}^2\right)}^{\left(n\right)}+\mathrm{De}{\overline{\lambda }}^{\left(n\right)}\frac{{\mathit{\tau }}^{(n+1)}-{\mathit{\tau }}^{\left(n\right)}}{\Delta t}=\mathrm{explicit}\mathrm{terms}.$$

This leads to the direct solution

$$(\mathrm{De}{\overline{\lambda }}^{\left(n\right)}+\xi \Delta t){\mathit{\tau }}^{(n+1)}=\mathrm{explicit}\mathrm{terms}.$$

(18)

where the respective $\mathrm{explicit}\mathrm{terms}$ are

$$\begin{array}{cc}\hfill [\mathrm{De}{\overline{\lambda }}^{\left(n\right)}-(1-\xi )\Delta t]{\tau }_{11}^{\left(n\right)}& +\Delta t\mathrm{De}{\overline{\lambda }}^{\left(n\right)}\left[{\tau }_{12}^{\left(n\right)}\frac{\partial }{\partial y}{u}^{\left(n\right)}+{\tau }_{11}^{\left(n\right)}\frac{\partial }{\partial y}log\left(1+\alpha T^{\left(n\right)}\right)\right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -\epsilon \Delta t({\tau }_{11}^2+{\tau }_{12}^2),\hfill \\ \hfill [\mathrm{De}{\overline{\lambda }}^{\left(n\right)}-(1-\xi )\Delta t]{\tau }_{12}^{\left(n\right)}& +\Delta t\mathrm{De}{\overline{\lambda }}^{\left(n\right)}\left[{\tau }_{22}^{\left(n\right)}\frac{\partial }{\partial y}{u}^{\left(n\right)}+{\tau }_{12}^{\left(n\right)}\frac{\partial }{\partial y}log\left(1+\alpha T^{\left(n\right)}\right)\right]\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & +\Delta t{\left({\eta }_p\right)}_{nf}\frac{\partial }{\partial y}{u}^{\left(n\right)}-\epsilon \Delta t({\tau }_{11}{\tau }_{12}+{\tau }_{12}{\tau }_{22}),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill [\mathrm{De}{\overline{\lambda }}^{\left(n\right)}-(1-\xi )\Delta t]{\tau }_{22}^{\left(n\right)}& +\Delta t\mathrm{De}{\overline{\lambda }}^{\left(n\right)}{\tau }_{22}^{\left(n\right)}\frac{\partial }{\partial y}log\left(1+\alpha T^{\left(n\right)}\right)-\epsilon \Delta t({\tau }_{12}^2+{\tau }_{22}^2).\hfill \end{array}$$

(21)

3.1. Graphical and Qualitative Results

Unless otherwise indicated, the following list of default values will be assumed for the various quantities and parameters:

$$\begin{array}{cc}\hfill & \alpha =0.01,\mathrm{Br}=1,\mathrm{Re}=1,\mathrm{Pr}=1,\mathrm{De}=2,\gamma =0.5,\beta =0.2,\Delta y=0.01,\hfill \\ \hfill & \Delta t=0.01,t=50,{\delta }_1=0.5,\zeta =1,\varphi =0.04,\epsilon =0.5,A_2=0.2,\chi =3,\xi =1.\hfill \end{array}$$

(22)

3.2. Time Development of Steady Smooth Solutions

Figure 2, Figure 3, Figure 4 and Figure 5 show the development (in time) of flow quantities until steady states are reached. Specifically, Figure 2 illustrates the transient development of nanofluid velocity, u, and nanofluid temperature, T. The fluid flow is driven by shear forces acting longitudinally on the upper wall. Initially, the bulk fluid is at rest. We notice from Figure 2—as expected—that the fluid velocity develops steadily in time from the initial state until the linear shear flow profile is achieved. Similarly, the fluid temperature is initially at zero throughout the flow geometry. As the flow develops in time, the exothermic reactions and mechanical dissipation processes taking place in the bulk flow lead to inevitable temperature increases in the bulk fluid. The channel walls are otherwise maintained at isothermal temperatures. We therefore notice from Figure 2—as expected—that the bulk flow temperature gradually increases in time until steady state parabolic temperature profiles are reached.

Figure 3 and Figure 4 show the time development of polymer stresses, ${\tau }_{11},{\tau }_{22},{\tau }_{12}$, and, hence, also of the second normal stress difference $N_1={\tau }_{11}-{\tau }_{22}$. As with the velocity and temperature profiles, we notice that the polymer stresses converge to fixed (reproducible) steady state solutions. The fixed (reproducibility) nature of the converged steady solutions is specifically illustrated in Figure 5, where the same converged values are obtained irrespective of the subsequent time of additional computation. For example, once the solutions converge to steady-state, at say, time $t=20$, then extending the computations to say, time $t=100$, will not change (qualitatively and quantitatively) the steady solutions.

3.3. Mesh-Size and Time-Step and Convergence

Our numerical algorithms use central differences to approximate the space derivatives. This means that our numerical schemes are second-order accurate in space. The accuracy in time is determined by the choice of $\xi$. Our computational results are all obtained with $\xi =1$. This means that our numerical algorithms (with $\xi =1$) are essentially first-order accurate in time. Taking, for example, $\xi =0.5$, as in the Crank–Nicolson scheme, would increase the accuracy in time to second-order. The choice $\xi =0.5$ would however require small time-steps ($\Delta t\ll 1$) to ensure that the computations converge in time. For this reason, we employ $\xi =1$; hence, we are free to use large time-steps, as illustrated in Figure 6. As shown in Figure 7, our numerical and computational algorithms also converge even with moderate values of the mesh-size, $\Delta y\approx 1$.

Figure 6 and Figure 7 therefore demonstrate—as required—that our numerical algorithms are largely independent of both mesh-size and time-step. Specifically, the computational algorithms efficiently reproduce the expected solutions for a broad range of mesh sizes and time-steps.

3.4. Development of Shear Banding

As noted in the introduction, shear banding phenomena in the shear flow of viscoelastic fluids represent observable and physical discontinuities in the shear rate profiles of the flow velocity. Shear banding phenomena have been observed in the shear flow of polymeric fluids of various types, including, for example, worm-like micelle solutions, polymer solutions, foams, telechelic polymers, granular flows, soft glasses, and polymer melts.

As noted in the introduction, there are two main scientifically-agreed mechanisms for shear banding, namely flow-induced inhomogeneities and constitutive instabilities. The Rolie-Poly viscoelastic constitutive model has been developed to explain the shear banding phenomena via flow-induced inhomogeneities. The Johnson-Segalman and Giesekus viscoelastic constitutive models allow for the mechanisms of the shear banding phenomena via constitutive instabilities if certain values of the viscoelastic material parameters are taken.

Specifically, when $\beta =0.95$ and $\epsilon =2$ in the Giesekus constitutive model, we note that the steady solutions (for the velocity u, temperature T, polymer-stresses ${\tau }_{11}$, ${\tau }_{12}$, ${\tau }_{22}$, and first normal stress difference $N_1$) exhibit shear banding phenomena, see Figure 8.

As indicated, and as in [27], the shear banding phenomena depicted in Figure 8 is linked to the viscoelastic parameters, specifically the polymer-to-total viscosity ratio ($\beta$) and the Giesekus nonlinear parameter ($\epsilon$) in the present study. In [27], the shear banding was linked to $\beta$ and the Johnson-Segalman parameter. Taking $\epsilon =0$ reduces the Giesekus viscoelastic constitutive model to the Oldroyd-B model, which is not capable of exhibiting shear banding, as illustrated in Figure 9 for the velocity u, temperature T, polymer-stresses ${\tau }_{11}$, ${\tau }_{12}$, ${\tau }_{22}$, and first normal stress difference $N_1$. Similarly, varying the values of $\beta$ shows that the onset and prevalence of shear banding is linked to higher values of $\beta$, as long as the value of $\epsilon$ is high enough, see Figure 10.

The phenomena of shear banding via constitutive instabilities is extensively investigated in [27] using the Johnson-Segalman viscoelastic constitutive model. In particular, it is conclusively demonstrated in [27] that shear banded solutions are essentially weak solutions of the governing partial differential equations for viscoelastic fluid flow and that an infinite number of such solutions exist for any set of shear banding material parameters. However, the results in [27] illustrate that, for each set of shear banding material parameter values, the computational results will always reproduce the same unique (albeit weak) solution! The unique selected shear rate path has been scientifically explained using the concept of stress diffusion as demonstrated in [28].

Figure 11 and Figure 12 illustrate the transient development of steady shear banded solutions, demonstrating, as in [27,28] the robustness and reproducibility of the shear banded solutions.

3.5. Thermal Runway

In industrial applications involving exothermic reactions, it is fundamental to ensure thermal stability and avoid thermal runaway phenomena. Thermal runaway can be demonstrated by increasing the values of the exothermic reaction parameter, ${\delta }_1$. Figure 13 illustrates the effects of ${\delta }_1$ on the steady-state VFBN velocity and temperature, specifically the drastic increase of temperature with increasing ${\delta }_1$.

Figure 14 gives a comparative illustration of the thermal runaway phenomena for Newtonian fluid-based and viscoelastic fluid-based nanofluids, i.e., NFBN and VFBN. The results for the Giesekus type VFBN computations are consistent with those in [23], which are based on the particle-free Oldroyd-B fluid computations. In particular, our results demonstrate that NFBN are more susceptible to thermal runaway than are the VFBN. This is consistent with the results in [23]. Essentially, the current results and those in [23] demonstrate that fluid viscoelasticity fundamentally acts to reduce large temperature build-ups in the exothermic shear flow setup and, hence, can be used to mitigate against rapid and finite-time temperature blow-up, i.e., thermal runaway. The present study will additionally demonstrate that shear banding phenomena can further enhance these favorable HTR characteristics of viscoelastic fluids, thus extending the scope of the investigations in [1,23].

4. Parameter Dependence of Solutions under Shear Banding Conditions

The dependence of thermal runaway on ${\delta }_1$ and that of shear banding on $\beta$ and $\epsilon$ have already been illustrated.

We now demonstrate, more generally, the sensitivity of the flow variables (solutions) to variations in these and other embedded parameters under shear banding conditions, i.e., $\beta =0.95$, $\epsilon =2$.

Figure 15 shows the dependence of the VFBN temperature field and thermal conductivity on nanoparticle volume fraction, $\varphi$. As expected, both the VFBN temperature and thermal conductivity increase with increasing $\varphi$. This is clearly expected, since the higher the concentration of heat-conducting nanoparticles (i.e., higher $\varphi$), the higher the corresponding Therm-C capabilities of the resultant nanofluid, and hence, the higher the temperatures in the nanofluid.

Figure 16 shows the dependence of the VFBN velocity and temperature on the Brinkman number, Br. Being directly proportional to the heat sources in the energy equation, increases in the Brinkman number, Br, would correspondingly increase the fluid temperature. Figure 16 also illustrates that the Brinkman number would not affect the shear banded nature of the velocity profile.

As with the thermal runaway analysis presented earlier (see Figure 13 and Figure 14) as well as the comparative literature, say [23,27], Figure 17 shows, as expected, that VFBNs are more resistant to temperature increases than NFBN.

Figure 17 additionally shows, as already expected, that the nanofluid temperature increases with increasing Br, but it decreases with increasing polymer viscosity, $\beta$. This further illustrates the viability of viscoelasticity (for which the polymer viscosity is a proxy) in dampening unwanted temperature increases and temperature blow-ups. As in [23], Figure 17 shows that the increase of temperature, with respect to Br, is linear for viscoelastic flow as opposed to the exponential increase observed with respect to ${\delta }_1$. This therefore means that increases in Br are not expected to lead to thermal runaway in shear flow of VFBN.

We conclude by demonstrating the expected increase of both the VFBN temperature and thermal conductivity with increasing activation energy parameter, $\alpha$, see Figure 18. As would be logically expected, the higher the activation energy strength, the faster the exothermic reaction processes and, hence, the higher the nanofluid thermal conductivities and obtainable temperatures, as illustrated in Figure 18.

5. Concluding Remarks

We employed a versatile and convergent numerical and computational algorithm, based on semi-implicit finite-difference methods to study the shear banding and thermal runaway characteristics of VFBN, which use Giesekus base fluids. We observed remarkable HTR and Therm-C enhancement in the VFBN under shear banding conditions of the material parameters. Specifically, we demonstrated that the VFBN are less susceptible to the thermal runaway than NFBN, and that this reduced susceptibility to thermal runaway phenomena is further enhanced under shear banding conditions—in particular, under increased polymer viscosity conditions.