# Stagnation Flow of a SWCNT Nanofluid towards a Plane Surface with Heterogeneous-Homogeneous Reactions

Authors: Sohail Ahmed; Hang Xu; Qiang Sun
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The homogeneous-heterogeneous reaction in the boundary layer flow of a water-based nanofluid in the stagnation-point region of a plane surface is investigated. The type of small particles explored here is the single-walled carbon nanotubes. The homogeneous nanofluid model is employed for description of behaviours of nanofluids. Here, the homogeneous (bulk) reaction is isothermal cubic autocatalytic, while the heterogeneous (surface) reaction is single, isothermal, and first order. The steady state of this system is analysed in detail, with equal diffusion coefficients being considered for both reactants and autocatalysts. Multiple solutions of the reduced system are captured for some particular sets of physical parameters, which seem to be overlooked in all previous published works with regard to studies of homogeneous-heterogeneous reactions modeled by homogeneous nanofluid models. Besides, we discover the significant limitation of previous conclusion about that the solutions by homogeneous nanofluid flow models can be recovered from those by regular fluids.


---

## Body

## 1. Introduction

Chemically reacting systems in liquids or biochemical systems such as catalysis and burning admit both homogeneous and heterogeneous reactions. However, the correlation between homogeneous and heterogeneous reactions is rather complex to handle. The reaction heat that is generated or absorbed during the chemical reactions has strong influence on the flow and heat transfer in the surrounding fluid that in turn affects the local concentrations of the reactants and the products and eventually the performance of the chemical reactions.Motivated by experimental observations of Williams et al. [1] and Song et al. [2] with regard to homogeneous-heterogeneous reactions on methane/ammonia and propane oxidation over a platinum, Chaudhary and Merkin [3, 4] developed a mathematical model on homogeneous-heterogeneous reactions in which the bulk reaction is expressed by the isothermal cubic autocatalytic kinetics and the surface reaction is governed by the first-order kinetics. Xu [5] made an extension to the reaction model of Chaudhary and Merkin [3, 4] to introduce the effect of the homogeneous-heterogeneous reactions into the energy equation.It has been known that heat transfer rate of base fluids could be enhanced significantly as highly conductive solid particles are added into them [6–8]. For example, Choi et al. [6] discovered the prominent heat transfer enhancement of working fluid when copper nanotubes are in suspension of water. As such, if the nanoparticles are added in the working fluid to make it become a nanofluid, since the thermophysical properties of the fluid have been changed, the chemical reaction happening in that fluid can be affected. When assuming that the suspended nanoparticles are not involved in the chemical reactions of interest, Xu and his collaborators [9, 10] further extended the homogeneous-heterogeneous reaction model to nanofluids by using Buongiorno’s approach [11]. Hayat et al. [12] applied the homogeneous-heterogeneous reaction model to three-dimensional flow of a nanofluid in the suspension of carbon nanotubes over a nonlinear stretchable surface. Bachok et al. [13] studied a stagnation region flow towards a stretching sheet with homogeneous-heterogeneous reaction effects. Other works engaged in the application of the homogeneous-heterogeneous model to different flow and heat transfer problems can be found in the literature studies [14–17].The suspended nanoparticles in nanofluids can be either metal, metal oxide, carbon nanotubes, or the combinations of them. Recently, carbon nanotubes (CNTs) have attracted much attention since they were discovered by Iijima [18] through discharging carbon electrodes to develop the new multiwalled carbon nanotubes (MWCNTs). Bethune et al. [19] then reported the development of single-walled carbon nanotubes (SWCNTs) that are structurally different from the multiwalled carbon nanotubes (MWCNTs). Since then, carbon nanotubes have become an attractive topic owing to their unique mechanical and thermal properties. Ramasubramaniam et al. [20] fabricated homogeneous single-wall carbon nanotubes (SWCNTs) for electrical conductivity applications. Xue [21] established a relation for effective thermal conductivity of carbon nanotubes based on composites. The high diffusion conductivity of the single-wall carbon nanotubes was then checked by some researchers such as Ma et al. [22], Khoshmehr et al. [23], and Park and Kim [24].As aforementioned, chemical reactions in liquids are complex, and different reaction rates can happen even under the same conditions. From the modelling point of view, this means that the solutions of the theoretical model for the chemical reaction of interest may not be unique. Though a few studies on homogeneous-heterogeneous chemical reactions in the boundary layer flow and heat transfer based on the homogeneous nanofluid model [8] have been done, no solution multiplicities yet have been reported. Based on the discussions by Chaudhary and Merkin [3, 4] as well as Xu and his collaborators [5, 9, 10], multiple solutions could be possible for homogeneous-heterogeneous chemical reactions with the homogeneous nanofluid model. In this paper, we investigated the homogeneous-heterogeneous reactions in a nanofluid in suspension of the single-walled carbon nanotubes (SWCNTs) when equal diffusion coefficients are considered for both autocatalyst and reactants. We found the existence of solution bifurcations which, to the best of our knowledge, has not been reported before. We also testified the validation of the previous conclusion [25] on that the solutions by homogeneous nanofluid flow models can be recovered from those for regular fluids and commented that it is only valid for a few special cases. This work can be deemed as a significant improvement of previous studies on homogeneous-heterogeneous reactions in a homogeneous nanofluid.


## 2. Mathematical Modeling

### 2.1. Problem Description

Consider a steady boundary layer flow and heat transfer of a SWCNT-nanofluid in the stagnation-point region towards a plane surface in the presence of homogeneous-heterogeneous chemical reactions. The homogeneous-heterogeneous reaction model proposed by Chaudhary and Merkin [3, 4] is applied to describe the behaviours of chemical reactions. In this model, the homogeneous reaction takes place in the bulk, which is isothermal cubic autocatalytic and governed by(1)A+2B⟶3B,rate=k1ab2,while the heterogeneous reaction occurs on the catalyst surface, which is single, isothermal, and the first order and written by(2)A⟶B,rate=ksa.In the aforementioned equations,a and b denote concentrations of the species A and B, respectively, and k1 and ks are the homogeneous and the heterogeneous reaction rates, respectively. Note that Williams et al. [1] and Song et al. [2] have confirmed the existence of homogeneous-heterogeneous reactions via experiments about methane/ammonia and propane oxidation over platinum.It is known that the homogeneous flow model [8] has been used to predict the behaviours of homogeneous-heterogeneous chemical reactions in nanofluids’ flow. By this approach, the thermal and physical properties of a nanofluid are derived from the base fluid and the suspended solid nanoparticles [21] as(3a)μnf=μf1−ϕ2.5,(3b)νnf=μfρnf,(3c)ρnf=1−ϕρf+ϕρCNT,(3d)ρcpnf=1−ϕρcpf+ϕρcpCNT,(3e)αnf=knfρcpnf,(3f)knf=1−ϕ+2ϕkCNT/kCNT−kflnkCNT+kf/2kf1−ϕ+2ϕkf/kCNT−kflnkCNT+kf/2kfkf,where ϕ is the volume fraction of nanoparticles, μ is the dynamic viscosity, ν is the kinematic viscosity, ρ is the density, α is the thermal diffusivity, cp is the heat capacity, and k is the thermal conductivity. The subscripts nf, f, and CNT denote nanofluid, fluid, and carbon nanotube (SWCNT), respectively. The thermophysical properties of water and SWCNTs used in this study are listed in Table 1.Table 1
Thermophysical characteristics of base fluid and SWCNTs.
Physical characteristics
Base fluid (water)
Nanoparticles (SWCNTS)
ρ
k
g
/
m
3
997.1
2600
c
p
J
/
k
g
k
4179
425
k
W
/
mk
0.613
6600Correspondingly, the governing equations describing the conservations of total mass, momentum, and energy, as well as chemical reaction diffusion, in the presence of single-wall carbon nanotubes (SWCNTs), in the framework of the boundary-layer approximations, are written as(4a)∂u∂x+∂v∂y=0,(4b)u∂u∂x+v∂u∂y=νnf∂2u∂y2+ue∂ue∂x,(4c)u∂T∂x+v∂T∂y=αnf∂2T∂y2,(4d)u∂a∂x+v∂a∂y=DA∂2a∂y2−k1ab2,(4e)u∂b∂x+v∂b∂y=DB∂2b∂y2+k1ab2,subject to the boundary conditions(5a)u=0,v=0,T=Tw,DA∂a∂y=−DB∂b∂y=ksa,aty=0,(5b)u=ue=cx,T⟶T∞,a⟶a∞,b⟶0,asy⟶∞,where u and v are the fluid velocity components in the x- and y-directions, respectively, ue is the incoming velocity of external flow, DA and DB are the diffusion coefficient of species A and B, respectively, Tw is the surface temperature, T∞ is the reference temperature, c is a constant, and a∞ is the concentration of species A at the far field.


### 2.2. Scaling Transformation

In the studies of flow and heat transfer in the boundary layers, it is a common practice to introduce the boundary layer assumptions [26] to reduce the governing partial differential equations to ordinary differential equations via similarity transformation. The resulting equations contain self-similar solutions which are similar to themselves as independent and dependent variables are appropriately scaled.We utilize the same idea to handle the problem considered in this work. In doing so, we introduce the following similarity variables:(6)η=cνf1/2y,u=cxf′η,v=−cνf1/2fη,θη=T−T∞Tw−T∞,ψη=aa∞,χη=ba∞.Substituting the similarity variables (6) into equation (4), continuity equation (4a) is automatically satisfied, and the rest equations are reduced to(7a)11−ϕ2.51−ϕ+ϕρCNT/ρff‴+ff″−f′2+1=0,(7b)11−ϕ+ϕρcpCNT/ρcpfknfkfθ″+Prfθ′=0,(7c)1ScAψ″+fψ′−Kψχ2=0,(7d)ϵScBχ″+fχ′+Kψχ2=0,subject to the boundary conditions(8a)f0=0,f′0=0,θ0=1,ψ′0=−ϵχ′0=Ksψ0,(8b)f′∞=1,θ∞=0,ψ∞=1,χ∞=0,where Pr=νf/αf is the Prandtl number, K=k1a∞2/c is the strength coefficient of the homogeneous reaction, ScA=νf/DA and ScB=νf/DB are the Schmidt numbers, ϵ=DB/DA is the ratio between the diffusion coefficients, and Ks=ksc/νf−1/2/DA is the strength coefficient due to the heterogeneous reaction. Note that, here, all those physical parameters are defined based on the properties of the base fluid so that the influence of nanoparticles can be clearly reflected.In many practical applications, the diffusion of speciesA and B is in comparable sizes, which indicates that the diffusion coefficients DA and DB can be deemed as equal [3]. As a result, we are able to choose ϵ=1. In this case, the following relation holds:(9)ψη+χη=1.Using equations (9), (7c), and (7d) which are combined to obtain(10)1ScAψ″+fψ′−Kψ1−ψ2=0.


### 2.3. Remarks on the Homogeneous Flow Model

It is worth mentioning to this end that Magyari [25] once argued that the solutions obtained by the homogeneous nanofluid flow model could be recovered by the corresponding regular fluid’s results, and no further efforts need to be paid for investigation of homogeneous modelling nanofluid flows. Here, we testify the validity of his conclusion via equations (7a), (7b), and (10). In doing so, we introduce the following scaling transformations:(11)fη=λfFξ,θη=λθΘξ,ψη=λψΨξ,η=ληξ.where λf, λθ, λψ, λχ, and λη are positive constants.Substituting equation (11) into equations (7a), (7b), and (10), we obtain(12a)λfλη31−ϕ2.51−ϕ+ϕρCNT/ρfF‴+λf2λη2FF″−λf2λη2F′2+1=0,(12b)λη1−ϕ+ϕρcpCNT/ρcpfknfkfΘ″+λfPrFΘ′=0,(12c)λη2ScAΨ″+λfληfΨ′−KΨ1−2λψΨ+λψ2Ψ2=0.To resume the governing equations to Newtonian fluid’s ones, it must hold(13a)λfλη31−ϕ2.51−ϕ+ϕρCNT/ρf=λf2λη2=1,(13b)λη1−ϕ+ϕρcpCNT/ρcpfknfkf=λf,(13c)λη2=λfλη=1=λψ=λψ2.From equation (13a), we obtain(14)λf=11−ϕ2.51−ϕ+ϕρCNT/ρf1/2,λfλη=1.However, from equation (13c), we have(15)λf=λη=λψ=1.If only the flow field is considered, then the fixed relation between homogeneous modelling nanofluid and regular fluid’s solutions is readily discovered from correlation (14). However, as multiple physical fields are taken into account, such fixed relation is no longer in existence. This is to say, it is impossible to recover homogeneous modelling nanofluid’s solutions via simple scaling transformations. Therefore, we are able to conclude that Magyari’s comment [25] is only available for a few specific problems. For most fluid flow problems, especially those involving into multiphysical fields, the homogeneous nanofluid model could not be replaced by scaling transformations so that it plays a significant role on investigation of flow behaviours of dilute nanofluids.


### 2.4. Computational Method

Solutions of governing equations (7a), (7b), and (10) subject to the boundary conditions (8) are solved by means of the shooting technique [27]. The integration distances are chosen from a range of values such as η∞=20,30,40 with the step size being set to Δη=0.005,0.01,0.05, respectively. If all computational results for a prescribed set of physical parameters with different integration distances and step sizes remain unchanged, then we think the solutions are stable. In this way, it is found that, for η∞=20 and Δη=0.01, the stability of multiple solutions can be kept. The Newtonian iteration is executed with the convergence criteria of 1.0×10−8 for all considered cases to examine the accuracy of results.


## 2.1. Problem Description

Consider a steady boundary layer flow and heat transfer of a SWCNT-nanofluid in the stagnation-point region towards a plane surface in the presence of homogeneous-heterogeneous chemical reactions. The homogeneous-heterogeneous reaction model proposed by Chaudhary and Merkin [3, 4] is applied to describe the behaviours of chemical reactions. In this model, the homogeneous reaction takes place in the bulk, which is isothermal cubic autocatalytic and governed by(1)A+2B⟶3B,rate=k1ab2,while the heterogeneous reaction occurs on the catalyst surface, which is single, isothermal, and the first order and written by(2)A⟶B,rate=ksa.In the aforementioned equations,a and b denote concentrations of the species A and B, respectively, and k1 and ks are the homogeneous and the heterogeneous reaction rates, respectively. Note that Williams et al. [1] and Song et al. [2] have confirmed the existence of homogeneous-heterogeneous reactions via experiments about methane/ammonia and propane oxidation over platinum.It is known that the homogeneous flow model [8] has been used to predict the behaviours of homogeneous-heterogeneous chemical reactions in nanofluids’ flow. By this approach, the thermal and physical properties of a nanofluid are derived from the base fluid and the suspended solid nanoparticles [21] as(3a)μnf=μf1−ϕ2.5,(3b)νnf=μfρnf,(3c)ρnf=1−ϕρf+ϕρCNT,(3d)ρcpnf=1−ϕρcpf+ϕρcpCNT,(3e)αnf=knfρcpnf,(3f)knf=1−ϕ+2ϕkCNT/kCNT−kflnkCNT+kf/2kf1−ϕ+2ϕkf/kCNT−kflnkCNT+kf/2kfkf,where ϕ is the volume fraction of nanoparticles, μ is the dynamic viscosity, ν is the kinematic viscosity, ρ is the density, α is the thermal diffusivity, cp is the heat capacity, and k is the thermal conductivity. The subscripts nf, f, and CNT denote nanofluid, fluid, and carbon nanotube (SWCNT), respectively. The thermophysical properties of water and SWCNTs used in this study are listed in Table 1.Table 1
Thermophysical characteristics of base fluid and SWCNTs.
Physical characteristics
Base fluid (water)
Nanoparticles (SWCNTS)
ρ
k
g
/
m
3
997.1
2600
c
p
J
/
k
g
k
4179
425
k
W
/
mk
0.613
6600Correspondingly, the governing equations describing the conservations of total mass, momentum, and energy, as well as chemical reaction diffusion, in the presence of single-wall carbon nanotubes (SWCNTs), in the framework of the boundary-layer approximations, are written as(4a)∂u∂x+∂v∂y=0,(4b)u∂u∂x+v∂u∂y=νnf∂2u∂y2+ue∂ue∂x,(4c)u∂T∂x+v∂T∂y=αnf∂2T∂y2,(4d)u∂a∂x+v∂a∂y=DA∂2a∂y2−k1ab2,(4e)u∂b∂x+v∂b∂y=DB∂2b∂y2+k1ab2,subject to the boundary conditions(5a)u=0,v=0,T=Tw,DA∂a∂y=−DB∂b∂y=ksa,aty=0,(5b)u=ue=cx,T⟶T∞,a⟶a∞,b⟶0,asy⟶∞,where u and v are the fluid velocity components in the x- and y-directions, respectively, ue is the incoming velocity of external flow, DA and DB are the diffusion coefficient of species A and B, respectively, Tw is the surface temperature, T∞ is the reference temperature, c is a constant, and a∞ is the concentration of species A at the far field.


## 2.2. Scaling Transformation

In the studies of flow and heat transfer in the boundary layers, it is a common practice to introduce the boundary layer assumptions [26] to reduce the governing partial differential equations to ordinary differential equations via similarity transformation. The resulting equations contain self-similar solutions which are similar to themselves as independent and dependent variables are appropriately scaled.We utilize the same idea to handle the problem considered in this work. In doing so, we introduce the following similarity variables:(6)η=cνf1/2y,u=cxf′η,v=−cνf1/2fη,θη=T−T∞Tw−T∞,ψη=aa∞,χη=ba∞.Substituting the similarity variables (6) into equation (4), continuity equation (4a) is automatically satisfied, and the rest equations are reduced to(7a)11−ϕ2.51−ϕ+ϕρCNT/ρff‴+ff″−f′2+1=0,(7b)11−ϕ+ϕρcpCNT/ρcpfknfkfθ″+Prfθ′=0,(7c)1ScAψ″+fψ′−Kψχ2=0,(7d)ϵScBχ″+fχ′+Kψχ2=0,subject to the boundary conditions(8a)f0=0,f′0=0,θ0=1,ψ′0=−ϵχ′0=Ksψ0,(8b)f′∞=1,θ∞=0,ψ∞=1,χ∞=0,where Pr=νf/αf is the Prandtl number, K=k1a∞2/c is the strength coefficient of the homogeneous reaction, ScA=νf/DA and ScB=νf/DB are the Schmidt numbers, ϵ=DB/DA is the ratio between the diffusion coefficients, and Ks=ksc/νf−1/2/DA is the strength coefficient due to the heterogeneous reaction. Note that, here, all those physical parameters are defined based on the properties of the base fluid so that the influence of nanoparticles can be clearly reflected.In many practical applications, the diffusion of speciesA and B is in comparable sizes, which indicates that the diffusion coefficients DA and DB can be deemed as equal [3]. As a result, we are able to choose ϵ=1. In this case, the following relation holds:(9)ψη+χη=1.Using equations (9), (7c), and (7d) which are combined to obtain(10)1ScAψ″+fψ′−Kψ1−ψ2=0.


## 2.3. Remarks on the Homogeneous Flow Model

It is worth mentioning to this end that Magyari [25] once argued that the solutions obtained by the homogeneous nanofluid flow model could be recovered by the corresponding regular fluid’s results, and no further efforts need to be paid for investigation of homogeneous modelling nanofluid flows. Here, we testify the validity of his conclusion via equations (7a), (7b), and (10). In doing so, we introduce the following scaling transformations:(11)fη=λfFξ,θη=λθΘξ,ψη=λψΨξ,η=ληξ.where λf, λθ, λψ, λχ, and λη are positive constants.Substituting equation (11) into equations (7a), (7b), and (10), we obtain(12a)λfλη31−ϕ2.51−ϕ+ϕρCNT/ρfF‴+λf2λη2FF″−λf2λη2F′2+1=0,(12b)λη1−ϕ+ϕρcpCNT/ρcpfknfkfΘ″+λfPrFΘ′=0,(12c)λη2ScAΨ″+λfληfΨ′−KΨ1−2λψΨ+λψ2Ψ2=0.To resume the governing equations to Newtonian fluid’s ones, it must hold(13a)λfλη31−ϕ2.51−ϕ+ϕρCNT/ρf=λf2λη2=1,(13b)λη1−ϕ+ϕρcpCNT/ρcpfknfkf=λf,(13c)λη2=λfλη=1=λψ=λψ2.From equation (13a), we obtain(14)λf=11−ϕ2.51−ϕ+ϕρCNT/ρf1/2,λfλη=1.However, from equation (13c), we have(15)λf=λη=λψ=1.If only the flow field is considered, then the fixed relation between homogeneous modelling nanofluid and regular fluid’s solutions is readily discovered from correlation (14). However, as multiple physical fields are taken into account, such fixed relation is no longer in existence. This is to say, it is impossible to recover homogeneous modelling nanofluid’s solutions via simple scaling transformations. Therefore, we are able to conclude that Magyari’s comment [25] is only available for a few specific problems. For most fluid flow problems, especially those involving into multiphysical fields, the homogeneous nanofluid model could not be replaced by scaling transformations so that it plays a significant role on investigation of flow behaviours of dilute nanofluids.


## 2.4. Computational Method

Solutions of governing equations (7a), (7b), and (10) subject to the boundary conditions (8) are solved by means of the shooting technique [27]. The integration distances are chosen from a range of values such as η∞=20,30,40 with the step size being set to Δη=0.005,0.01,0.05, respectively. If all computational results for a prescribed set of physical parameters with different integration distances and step sizes remain unchanged, then we think the solutions are stable. In this way, it is found that, for η∞=20 and Δη=0.01, the stability of multiple solutions can be kept. The Newtonian iteration is executed with the convergence criteria of 1.0×10−8 for all considered cases to examine the accuracy of results.


## 3. Solution Dependence onK

It has been already known from Chaudhary and Merkin [3] that the strength coefficients of homogeneous and heterogeneous reactions are key factors to generate solution bifurcations. We therefore testify the roles of K and Ks on the solution behaviours with various nanoparticle volume fractions. In our computation, the Prandtl number is fixed, Pr=6.8, for the base fluid (water) throughout the study. Since the homogeneous nanofluid model is only suitable for the dilute nanofluid, the range of the volume fraction of nanoparticles is chosen from 0 to 0.2 0≤ϕ≤0.2. Particularly, from the definitions of nanofluids’ parameters denoted in equation (3), it is readily known that the case ϕ=0 corresponds to the regular fluid.


### 3.1. Solutions for Small and ModerateK

We start our discussion by considering the variation ofψ0 due to the jointly effects of K and Ks in the presence of SWCNT nanoparticles. Clearly, the multiple solutions still exist as the nanoparticles are considered. For example, for ϕ=0.1 and ScA=1, multiple solutions are found for Ks=0.025 and Ks=0.05 as K in particular regions. However, for Ks=0.1, only unique solution can be found, as shown in Figure 1. It is also observed from this figure that the hysteresis point is obtained for the decreasing values of Ks (say KSH). For all values, Ks<KSH, multiple solution are possible.Figure 1
ψ
0 against K for some values of Ks with Pr=6.8, ϕ=0.1, and ScA=1.0.The bifurcation points for possible ranges ofK are listed in Table 2. It is seen in this table that multiple solutions are found in the ranges Kmin≤K≤Kmax for Ks=0.025 and Kmin≤K≤Kmax for Ks=0.05, while only one unique solution is discovered for Ks=0.1.Table 2
Minimum and maximum values (sayKmin and Kmax) of K at several prescribed values of Ks in the case of ϕ=0.1, ScA=1.0, and Pr=6.8.
Bifurcation point
K
s
=
0.01
K
s
=
0.025
K
s
=
0.1
K
min
2.696
2.659
—
K
max
14.922
6.505
—We then consider the effect of the nanoparticle volume fractionϕ on ψ0 for different values of K. In doing so, we choose Ks=0.025, Pr=6.8, and ScA=1 as an illustrative case where multiple solutions are available. As shown in Figure 2, as ϕ increases, the bifurcation points for both maximum and minimum shift towards the left side. This indicates that the addition of nanoparticles is helpful to keep the stability of chemical reaction for large K, while it could lead to the instability of chemical reaction for small K. The maximum and minimum points for solution multiplicities are listed in Table 3. It is clearly shown that multiple solution regions are 2.723≤K≤6.738 for ϕ=0, 2.659≤K≤6.505 for ϕ=0.1, and 2.568≤K≤6.174 for ϕ=0.2.Figure 2
ψ
0 against K for some values of ϕ with Pr=6.8, Ks=0.025, and ScA=1.0.Table 3
Minimum and maximum values ofK at several prescribed values of ϕ in the case of Ks=0.025, ScA=1.0, and Pr=6.8.
Bifurcation point
ϕ
=
0
ϕ
=
0.1
ϕ
=
0.2
K
min
2.723
2.659
2.568
K
max
6.738
6.505
6.174The influence ofKs on ψ0 for different values of K at ϕ=0.1, ScA=1, and Pr=6.8 is presented in Figure 3. It is observed that only one unique solution exists as the homogeneous reaction rate K=2, while multiple solutions are captured for K=4 and K=8 provided that Ks is properly chosen. It is worth mentioning that, in our model assumption, the heterogeneous reaction rate Ks is expected to be always positive. In the meantime, ψη has to be equal or greater than zero and less than or equal to 1 (owing to a≤a∞). This means that all solutions for Ks<0 and ψ>1 (or ψ<0) are physically meaningless. Particularly, the system can keep the self-sustained surface reaction without the homogeneous reaction being involved K=0 when nanoparticles are appended. As the hysteresis point is reached for increasing value of K (say KH), multiple solutions are possible for all values of K>KH. Figure 4 shows the effect of Ks on ψ0 for some values of nanoparticle volume fraction ϕ. From the figure, it is observed that the variation of ϕ could not alter the behavior of solution multiplicity. In other words, the existence of multiple solutions is not affected by changes of nanoparticle volume fraction. The role of the nanoparticle volume fraction is that it can accelerate or postpone the chemical reaction rate.Figure 3
ψ
0 against Ks for some values of K with Pr=6.8, ϕ=0.1, and ScA=1.0.Figure 4
ψ
0 against Ks for some values of ϕ with Pr=6.8, K=1, and ScA=1.0.


### 3.2. Solutions for LargeK

It has been reported by Chaudhary and Merkin [3] that there is a shear-layer-type structure developed by the concentration profile ψη as K continuously increases. The similar structure is found to be still in existence when nanofluids are considered. As shown in Figure 5(a), the reaction gradually compacts to a relatively thin region that increasingly departs from the surface as ϕ evolves. This figure also shows that, as ϕ enlarges, the thickness of the fully reacted-state region grows larger, which leads to the thickness of the unreacted-state region which becomes thinner.Figure 5
(a) Concentration profileψη plotted against η for increasing values of K when ScA=1, Pr=6.8, Ks=0.1, ϕ=0.0, ϕ=0.1, and ϕ=0,2; (b) a plot of R(K), the position of the reaction front, for K larger obtained from equation (27) for ϕ=0.0, ϕ=0.1, and ϕ=0.2.
(a)
(b)From equation (10), it is expected that the asymptotic behavior for large K is Ks in the order of O(1). In the reaction zone, this equation suggests(16)η=RK+K−1/2Y,where RK is a constant and Y denotes large η. Here, the reaction front has a thickness of OK−1/2 and is centered on η=RK≫1. Therefore, RK is expanded as(17)RK=R0K−1/2+R1+⋯,asK⟶∞,where R0, R1, ⋯, are constants to be determined.Using equation (16) in equation (10), the following equation is obtained:(18)1ScAψ″+R+K−1/2Y−Δψ′−Kψ1−ψ2=0,subject to the boundary condition(19)ψ⟶1,Y⟶∞.Here, the prime denotes the differentiation with respect toY.Equation (18) is expected to have a solution of the form(20)ψY;K=ψ0Y+K−1/2ψ1Y+⋯.Substituting equation (20) into equation (18) and then omitting the high-order terms, at the first-order, it reads(21)1ScAψ0″+R0ψ0′−ψ01−ψ02=0.Taking boundary condition (19) into account, it is expected that ψ0 is exponentially small in K in the inner wall region such that(22)ψ0⟶0,asY⟶−∞.As briefly discussed in [28], equation (21) subject to boundary condition equation (22) is the cubic-Fisher wave equation that is essentially an eigenvalue problem for R0. Its solution has exponential decaying property at infinity as expected. The reason is that the solutions to equation (21) have algebraic properties of OY∗ as Y⟶∞, where Y∗=YScA1/2, which in turn bring about the perturbation of ORe−1/2 to ψ in the external flow.Using the method proposed by Billingham and Nadeem [28], the unique solution for equation (21) is given as(23)R0=12ScA,ψ0=1+exp−Y∗2−1.AtOK−1/2, it holds(24)ψ1+12ψ1′−1−4ψ0+3ψ02ψ1=Δ−R1ScAψ0′,subject to the boundary conditions(25)ψ1⟶0,asY⟶∞,where prime denotes the differentiation with respect to Y∗.Since equation (24) includes a complementary function ψ0′ which satisfies the required boundary conditions, an additional boundary condition has to be introduced to give a solution to the full equation. To meet this requirement, the following correlation is introduced (refer to [3]):(26)ScAΔ−R1∫−∞∞ψ0′2expY∗2dY∗=0,which gives R1=Δ. Therefore, RK can be expressed as(27)RK=K1/22ScA+Δ+⋯,asK⟶∞.However, as concluded by Chaudhary and Merkin [3], equation (27) is not a complete solution since the condition at η=0 does not match to that in equation (23). Therefore, they suggested to consider the condition at the inner edge of the reaction front where the inner solution ψ0 decays exponentially small in K. As such, the following relations hold:(28)ζ=K−1Y,ψ=expKMζ;K,where ζ<0 is in the order of O(1) for the region. Due to equation (28), it is readily known that ψ0 possesses exponentially decaying behaviors for large K. Therefore, we have(29)ψ0∼exp−M0K,asK⟶∞,M0=0.596.Note that this relation still holds when the effect of nanofluids is considered.It is worth mentioning to this end that the nanofluid is considered in our case; therefore, the nanoparticle volume fraction plays a role onRK. From equation (27), we notice that RK is closely related to Δ which is equal to 0.647, 0.685, and 0.744 for ϕ=0.0, ϕ=0.1, and ϕ=0.2, respectively. This indicates that RK enhances as ϕ evolves for prescribed values of K, as shown in Figure 5(b).


## 3.1. Solutions for Small and ModerateK

We start our discussion by considering the variation ofψ0 due to the jointly effects of K and Ks in the presence of SWCNT nanoparticles. Clearly, the multiple solutions still exist as the nanoparticles are considered. For example, for ϕ=0.1 and ScA=1, multiple solutions are found for Ks=0.025 and Ks=0.05 as K in particular regions. However, for Ks=0.1, only unique solution can be found, as shown in Figure 1. It is also observed from this figure that the hysteresis point is obtained for the decreasing values of Ks (say KSH). For all values, Ks<KSH, multiple solution are possible.Figure 1
ψ
0 against K for some values of Ks with Pr=6.8, ϕ=0.1, and ScA=1.0.The bifurcation points for possible ranges ofK are listed in Table 2. It is seen in this table that multiple solutions are found in the ranges Kmin≤K≤Kmax for Ks=0.025 and Kmin≤K≤Kmax for Ks=0.05, while only one unique solution is discovered for Ks=0.1.Table 2
Minimum and maximum values (sayKmin and Kmax) of K at several prescribed values of Ks in the case of ϕ=0.1, ScA=1.0, and Pr=6.8.
Bifurcation point
K
s
=
0.01
K
s
=
0.025
K
s
=
0.1
K
min
2.696
2.659
—
K
max
14.922
6.505
—We then consider the effect of the nanoparticle volume fractionϕ on ψ0 for different values of K. In doing so, we choose Ks=0.025, Pr=6.8, and ScA=1 as an illustrative case where multiple solutions are available. As shown in Figure 2, as ϕ increases, the bifurcation points for both maximum and minimum shift towards the left side. This indicates that the addition of nanoparticles is helpful to keep the stability of chemical reaction for large K, while it could lead to the instability of chemical reaction for small K. The maximum and minimum points for solution multiplicities are listed in Table 3. It is clearly shown that multiple solution regions are 2.723≤K≤6.738 for ϕ=0, 2.659≤K≤6.505 for ϕ=0.1, and 2.568≤K≤6.174 for ϕ=0.2.Figure 2
ψ
0 against K for some values of ϕ with Pr=6.8, Ks=0.025, and ScA=1.0.Table 3
Minimum and maximum values ofK at several prescribed values of ϕ in the case of Ks=0.025, ScA=1.0, and Pr=6.8.
Bifurcation point
ϕ
=
0
ϕ
=
0.1
ϕ
=
0.2
K
min
2.723
2.659
2.568
K
max
6.738
6.505
6.174The influence ofKs on ψ0 for different values of K at ϕ=0.1, ScA=1, and Pr=6.8 is presented in Figure 3. It is observed that only one unique solution exists as the homogeneous reaction rate K=2, while multiple solutions are captured for K=4 and K=8 provided that Ks is properly chosen. It is worth mentioning that, in our model assumption, the heterogeneous reaction rate Ks is expected to be always positive. In the meantime, ψη has to be equal or greater than zero and less than or equal to 1 (owing to a≤a∞). This means that all solutions for Ks<0 and ψ>1 (or ψ<0) are physically meaningless. Particularly, the system can keep the self-sustained surface reaction without the homogeneous reaction being involved K=0 when nanoparticles are appended. As the hysteresis point is reached for increasing value of K (say KH), multiple solutions are possible for all values of K>KH. Figure 4 shows the effect of Ks on ψ0 for some values of nanoparticle volume fraction ϕ. From the figure, it is observed that the variation of ϕ could not alter the behavior of solution multiplicity. In other words, the existence of multiple solutions is not affected by changes of nanoparticle volume fraction. The role of the nanoparticle volume fraction is that it can accelerate or postpone the chemical reaction rate.Figure 3
ψ
0 against Ks for some values of K with Pr=6.8, ϕ=0.1, and ScA=1.0.Figure 4
ψ
0 against Ks for some values of ϕ with Pr=6.8, K=1, and ScA=1.0.


## 3.2. Solutions for LargeK

It has been reported by Chaudhary and Merkin [3] that there is a shear-layer-type structure developed by the concentration profile ψη as K continuously increases. The similar structure is found to be still in existence when nanofluids are considered. As shown in Figure 5(a), the reaction gradually compacts to a relatively thin region that increasingly departs from the surface as ϕ evolves. This figure also shows that, as ϕ enlarges, the thickness of the fully reacted-state region grows larger, which leads to the thickness of the unreacted-state region which becomes thinner.Figure 5
(a) Concentration profileψη plotted against η for increasing values of K when ScA=1, Pr=6.8, Ks=0.1, ϕ=0.0, ϕ=0.1, and ϕ=0,2; (b) a plot of R(K), the position of the reaction front, for K larger obtained from equation (27) for ϕ=0.0, ϕ=0.1, and ϕ=0.2.
(a)
(b)From equation (10), it is expected that the asymptotic behavior for large K is Ks in the order of O(1). In the reaction zone, this equation suggests(16)η=RK+K−1/2Y,where RK is a constant and Y denotes large η. Here, the reaction front has a thickness of OK−1/2 and is centered on η=RK≫1. Therefore, RK is expanded as(17)RK=R0K−1/2+R1+⋯,asK⟶∞,where R0, R1, ⋯, are constants to be determined.Using equation (16) in equation (10), the following equation is obtained:(18)1ScAψ″+R+K−1/2Y−Δψ′−Kψ1−ψ2=0,subject to the boundary condition(19)ψ⟶1,Y⟶∞.Here, the prime denotes the differentiation with respect toY.Equation (18) is expected to have a solution of the form(20)ψY;K=ψ0Y+K−1/2ψ1Y+⋯.Substituting equation (20) into equation (18) and then omitting the high-order terms, at the first-order, it reads(21)1ScAψ0″+R0ψ0′−ψ01−ψ02=0.Taking boundary condition (19) into account, it is expected that ψ0 is exponentially small in K in the inner wall region such that(22)ψ0⟶0,asY⟶−∞.As briefly discussed in [28], equation (21) subject to boundary condition equation (22) is the cubic-Fisher wave equation that is essentially an eigenvalue problem for R0. Its solution has exponential decaying property at infinity as expected. The reason is that the solutions to equation (21) have algebraic properties of OY∗ as Y⟶∞, where Y∗=YScA1/2, which in turn bring about the perturbation of ORe−1/2 to ψ in the external flow.Using the method proposed by Billingham and Nadeem [28], the unique solution for equation (21) is given as(23)R0=12ScA,ψ0=1+exp−Y∗2−1.AtOK−1/2, it holds(24)ψ1+12ψ1′−1−4ψ0+3ψ02ψ1=Δ−R1ScAψ0′,subject to the boundary conditions(25)ψ1⟶0,asY⟶∞,where prime denotes the differentiation with respect to Y∗.Since equation (24) includes a complementary function ψ0′ which satisfies the required boundary conditions, an additional boundary condition has to be introduced to give a solution to the full equation. To meet this requirement, the following correlation is introduced (refer to [3]):(26)ScAΔ−R1∫−∞∞ψ0′2expY∗2dY∗=0,which gives R1=Δ. Therefore, RK can be expressed as(27)RK=K1/22ScA+Δ+⋯,asK⟶∞.However, as concluded by Chaudhary and Merkin [3], equation (27) is not a complete solution since the condition at η=0 does not match to that in equation (23). Therefore, they suggested to consider the condition at the inner edge of the reaction front where the inner solution ψ0 decays exponentially small in K. As such, the following relations hold:(28)ζ=K−1Y,ψ=expKMζ;K,where ζ<0 is in the order of O(1) for the region. Due to equation (28), it is readily known that ψ0 possesses exponentially decaying behaviors for large K. Therefore, we have(29)ψ0∼exp−M0K,asK⟶∞,M0=0.596.Note that this relation still holds when the effect of nanofluids is considered.It is worth mentioning to this end that the nanofluid is considered in our case; therefore, the nanoparticle volume fraction plays a role onRK. From equation (27), we notice that RK is closely related to Δ which is equal to 0.647, 0.685, and 0.744 for ϕ=0.0, ϕ=0.1, and ϕ=0.2, respectively. This indicates that RK enhances as ϕ evolves for prescribed values of K, as shown in Figure 5(b).


## 4. Autoignition

Autoignition takes place atKs=0. To decide the occurrence time of the autoignition, equation (10) subject to the following boundary conditions needs to be taken into account:(30)ψ′0=0,ψ⟶1,asη⟶∞.As shown in Figure6(a), the nanoparticle volume fraction ϕ has perceptible effects on ψ0. For each ϕ, there is a critical value of Ka. For K>Ka, two different branches are available when other parameters are prescribed. These solutions correspond to the two points in Figure 3, where Ks=0 for the case ϕ=0.1, ScA=1, and Pr=6.8. Though not shown in the figures, it is expected that the lower and the upper branches for ψ0 approach to 0 and 1, respectively, as K approaches to infinity.Figure 6
Autoignition: (a) a graph ofψ0 against K obtained by solving equation (10) with boundary condition equation (30) for ScA=1, Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, (b) a graph of Ka plotted against ScA for Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, and (c) a graph of ψ0 when Ka=K plotted against ScA.
(a)
(b)
(c)Particularly, forK⟶∞, equation (10) suggests that the upper branch has the following form:(31)ψ=1−K−1ψ∗,where ψ∗ is a scaling transformation function. If we let η be unscaled, then the following leading order problem is derived:(32)ψ∗″+fψ∗′+ψ∗2=0,subject to the boundary conditions(33)ψ∗′0=0,ψ∗⟶1,asη⟶∞.By solving this equation numerically, we obtainψ∗0=0.833 for ϕ=0.0, ψ∗0=0.825 for ϕ=0.1, and ψ∗0=0.818 for ϕ=0.2 when ScA=1 and Pr=6.8, respectively. As a result, we obtain the asymptotic expressions of ψ0 for K⟶∞ as(34)ψ0∼1−0.833K−1+⋯,ϕ=0.0,(35)ψ0∼1−0.825K−1+⋯,ϕ=0.1,(36)ψ0∼1−0.818K−1+⋯,ϕ=0.2.For example, asymptotic expression (35) is shown in Figure 6(a) by the dashed curve which is in good agreement with the numerical solution, especially for large K. Also, from our previous analysis, the behavior of the lower solution branch becomes exponentially small as K grows considerably large.As expected, the critical valueKa is dependent on ScA and ϕ as well. The dependence of Ka on ScA is shown in Figure 6(b), with the corresponding values of ψ0 being shown in Figure 6(c). From these figures, we notice that Ka reduces monotonously as ScA evolves. On the contrary, for prescribed values of ScA, Ka decreases gradually as ScA enlarges. Some critical values Ka for various values of ScA are listed in Table 4 in which the comparison with the results of Chaudhary and Merkin [3] is presented, and good agreement is found.Table 4
K
a against ScA for various values of ϕ at Pr=6.8.
Sc
A
K
a
Chaudhary and Merkin [3]
ϕ
=
0.0
ϕ
=
0.1
ϕ
=
0.2
ϕ
=
0.0
0.3
3.413
3.358
3.277
0.5
3.1589
3.0989
3.0108
1
2.783
2.719
2.627
2.785
1.5
2.556
2.492
2.401
2.5
2.273
2.211
2.123


## 5. Solution Dependance onScA

It is obvious thatScA can alter solution behaviours. We, therefore, discuss its influence on ψη when ScA is small and large, respectively.


### 5.1. Solution for LargeScA

To check the significance ofScA, we take K and Ks are both in the order of O(1). From equation (10), it is known that solutions could be possible for large ScA, provided that(37)ψη=1−ScA−1/3Bξ,ξ=ScA1/3η,where Bξ and ξ are the scaling function and variable, respectively.Equation (37) clearly shows that ψ0⟶1 as ScA approaches to infinity. The reaction region therefore has the thickness of OScA−1/3 for considerably large ScA. Substituting equation (37) in (10), it yields(38)B″+C02ξ2−ScA−1/3ξ36+⋯B′+KB21−ScA−1/3B=0,subject to the boundary conditions(39)B0=−Ks1−ScA−1/3B0,B∞=0.B can be expanded, based on equations (38) and (39), as(40)Bξ,ScA=B0ξ+ScA−1/3B1ξ.From equation (37), we further obtain(41)ψ0=1−ScA−1/3B00+⋯,asScA⟶∞,which indicates that there is no possibility for the existence of multiple solutions at large ScA when K and KS are both in the order of O(1).However, equation (41) will lose its valid range when Ks and K are not small. To handle this situation, we write(42)Ks=ScA1/3Ks∗,K=ScA−1/3K∗,where Ks∗ and K∗ both are O1 and ξ=ScAη. Equation (10) is then reduced to the following leading order equation:(43)ψ″+C02ξ2ψ′−K∗ψ1−ψ2=0,subject to the boundary conditions(44)ψ′0=Ks∗ψ0,ψ∞=1.The hysteresis bifurcation point can then be calculated by solving equation (43). For example, we obtain Ks∗=0.109 and K∗=3.220 at ϕ=0.0, Ks∗=0.064 and K∗=2.568 at ϕ=0.1, and Ks∗=0.039 and K∗=1.989 at ϕ=0.2 which further leads to(45a)KsH∼0.109ScA1/3,KH∼3.220ScA−1/3,ScA⟶∞,whenϕ=0.0.(45b)KsH∼0.064ScA1/3,KH∼2.568ScA−1/3,ScA⟶∞,whenϕ=0.1.(45c)KsH∼0.039ScA1/3,KH∼1.989ScA−1/3,ScA⟶∞,whenϕ=0.2.As shown in Figure7, the increase in ScA causes the decrease in KH (Figure 7(a)) but the increase in KsH (Figure 7(b)).Figure 7
Graphs of (a)KH plotted against Sc for Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, (b) KSH plotted against Sc for Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, and (c) KH plotted against KSH, the region where multiple steady-state solutions are possible are also indicated.
(a)
(b)
(c)


### 5.2. Solution for SmallScA

In this case, the reaction region is thicker than the region of the boundary layer flow. Since the thick outer region is of extentOScA−1/2, we introduce new variable ϱ=ScA1/2η and leave ψ0 unscaled, and then write(46)ψϱ,ScA=ψ0ϱ+ScA1/2ψ1ϱ+⋯.Substituting equation (46) into (10) and using the new scalings, at the leading order, we obtain(47)ψ0″+ϱψ0′−Kψ01−ψ02=0,ψ0∞=1.In the inner boundary layer region, we expandψ0ϱ at the leading order, obtaining(48)ψ0ϱ=D0+D1ϱ+D2ϱ2+⋯,where Di are constants which can be determined successively. Particularly, using equation (47), we obtain D2=1/2KD01−D02.At the order ofOScA1/2, ψ1ϱ is linear and satisfies ψ1∞=0; we therefore expand it, as ϱ approaches to 0, by(49)ψ1=E0+E1ϱ+⋯,where Ei are constants.On the contrary, based on equations (48) and (49), it suggests(50)ψη,ScA=ψ0η+ScA1/2ψ1η+⋯.Note that the inner region has to satisfy the boundary conditions (8b) at η=0. Therefore, the leading order equation gives(51)ψ0″=0,subject to the boundary conditions(52)ψ0′=Ksψ00,ψ0⟶D0,asη⟶∞.The solution of this equation isD0=0. At the order of OScA1/2, it holds(53)ψ1″=0,subject to the boundary conditions(54)ψ1′0=Ksψ10,ψ1∼D1η+E1,asη⟶∞,which give E1=D1/Ks. As a result, we obtain(55)ψ1=D1η+1Ks.Thus, the appropriate inner condition for equation (47) is(56)ψ00=0.The solutions of equation (47) subject to the inner condition (56) determine that D1=ψ0′0 from which we obtain(57)ψ0∼D1KsScA1/2+…,asScA⟶0.Equation (57) indicates that there is a singularity at Ks=0. To determine the nature of the hysteresis bifurcation for small ScA, we put(58)Ks=Ks∗ScA1/2.The solution at the leading order is simplified as(59)ψ0η=D0.However, at the order ofOScA1/2, equation (53) still holds and satisfies(60)ψ1′0=Ksψ00,ψ1∼D1η+E1,asη⟶∞,which give D1=Ks∗D0. So, we have(61)ψ1=Ks∗D0η+E1.To balance the outer region with the inner boundary condition, we obtain(62)ψϱ,ScA∼D0+Ks∗D0ϱ+ScA1/2E1+⋯,asϱ⟶0.We, therefore, are able to find the solution of equation (46) by solving leading order equation (47) subject to the inner boundary condition(63)ψ0′=Ks∗ψ0,onϱ⟶0.It is now expected that equation (47) subject to the inner boundary conditions (63) has multiple solutions. The hysteresis bifurcated point can be determined in a similar way. We obtain the critical values Ks∗H=0.129 and KH=4.712 for ϕ=0.0, Ks∗H=0.102 and KH=4.311 for ϕ=0.1, and Ks∗H=0.080 and KH=3.881 for ϕ=0.2, which further give(64a)KH⟶4.712,KsH∼0.129ScA1/2,asScA⟶0,whenϕ=0.0,(64b)KH⟶4.311,KsH∼0.102ScA1/2,asScA⟶0,whenϕ=0.1,(64c)KH⟶3.881,KsH∼0.080ScA1/2,asScA⟶0,whenϕ=0.2.


## 5.1. Solution for LargeScA

To check the significance ofScA, we take K and Ks are both in the order of O(1). From equation (10), it is known that solutions could be possible for large ScA, provided that(37)ψη=1−ScA−1/3Bξ,ξ=ScA1/3η,where Bξ and ξ are the scaling function and variable, respectively.Equation (37) clearly shows that ψ0⟶1 as ScA approaches to infinity. The reaction region therefore has the thickness of OScA−1/3 for considerably large ScA. Substituting equation (37) in (10), it yields(38)B″+C02ξ2−ScA−1/3ξ36+⋯B′+KB21−ScA−1/3B=0,subject to the boundary conditions(39)B0=−Ks1−ScA−1/3B0,B∞=0.B can be expanded, based on equations (38) and (39), as(40)Bξ,ScA=B0ξ+ScA−1/3B1ξ.From equation (37), we further obtain(41)ψ0=1−ScA−1/3B00+⋯,asScA⟶∞,which indicates that there is no possibility for the existence of multiple solutions at large ScA when K and KS are both in the order of O(1).However, equation (41) will lose its valid range when Ks and K are not small. To handle this situation, we write(42)Ks=ScA1/3Ks∗,K=ScA−1/3K∗,where Ks∗ and K∗ both are O1 and ξ=ScAη. Equation (10) is then reduced to the following leading order equation:(43)ψ″+C02ξ2ψ′−K∗ψ1−ψ2=0,subject to the boundary conditions(44)ψ′0=Ks∗ψ0,ψ∞=1.The hysteresis bifurcation point can then be calculated by solving equation (43). For example, we obtain Ks∗=0.109 and K∗=3.220 at ϕ=0.0, Ks∗=0.064 and K∗=2.568 at ϕ=0.1, and Ks∗=0.039 and K∗=1.989 at ϕ=0.2 which further leads to(45a)KsH∼0.109ScA1/3,KH∼3.220ScA−1/3,ScA⟶∞,whenϕ=0.0.(45b)KsH∼0.064ScA1/3,KH∼2.568ScA−1/3,ScA⟶∞,whenϕ=0.1.(45c)KsH∼0.039ScA1/3,KH∼1.989ScA−1/3,ScA⟶∞,whenϕ=0.2.As shown in Figure7, the increase in ScA causes the decrease in KH (Figure 7(a)) but the increase in KsH (Figure 7(b)).Figure 7
Graphs of (a)KH plotted against Sc for Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, (b) KSH plotted against Sc for Pr=6.8, ϕ=0.0, ϕ=0.1, and ϕ=0.2, and (c) KH plotted against KSH, the region where multiple steady-state solutions are possible are also indicated.
(a)
(b)
(c)


## 5.2. Solution for SmallScA

In this case, the reaction region is thicker than the region of the boundary layer flow. Since the thick outer region is of extentOScA−1/2, we introduce new variable ϱ=ScA1/2η and leave ψ0 unscaled, and then write(46)ψϱ,ScA=ψ0ϱ+ScA1/2ψ1ϱ+⋯.Substituting equation (46) into (10) and using the new scalings, at the leading order, we obtain(47)ψ0″+ϱψ0′−Kψ01−ψ02=0,ψ0∞=1.In the inner boundary layer region, we expandψ0ϱ at the leading order, obtaining(48)ψ0ϱ=D0+D1ϱ+D2ϱ2+⋯,where Di are constants which can be determined successively. Particularly, using equation (47), we obtain D2=1/2KD01−D02.At the order ofOScA1/2, ψ1ϱ is linear and satisfies ψ1∞=0; we therefore expand it, as ϱ approaches to 0, by(49)ψ1=E0+E1ϱ+⋯,where Ei are constants.On the contrary, based on equations (48) and (49), it suggests(50)ψη,ScA=ψ0η+ScA1/2ψ1η+⋯.Note that the inner region has to satisfy the boundary conditions (8b) at η=0. Therefore, the leading order equation gives(51)ψ0″=0,subject to the boundary conditions(52)ψ0′=Ksψ00,ψ0⟶D0,asη⟶∞.The solution of this equation isD0=0. At the order of OScA1/2, it holds(53)ψ1″=0,subject to the boundary conditions(54)ψ1′0=Ksψ10,ψ1∼D1η+E1,asη⟶∞,which give E1=D1/Ks. As a result, we obtain(55)ψ1=D1η+1Ks.Thus, the appropriate inner condition for equation (47) is(56)ψ00=0.The solutions of equation (47) subject to the inner condition (56) determine that D1=ψ0′0 from which we obtain(57)ψ0∼D1KsScA1/2+…,asScA⟶0.Equation (57) indicates that there is a singularity at Ks=0. To determine the nature of the hysteresis bifurcation for small ScA, we put(58)Ks=Ks∗ScA1/2.The solution at the leading order is simplified as(59)ψ0η=D0.However, at the order ofOScA1/2, equation (53) still holds and satisfies(60)ψ1′0=Ksψ00,ψ1∼D1η+E1,asη⟶∞,which give D1=Ks∗D0. So, we have(61)ψ1=Ks∗D0η+E1.To balance the outer region with the inner boundary condition, we obtain(62)ψϱ,ScA∼D0+Ks∗D0ϱ+ScA1/2E1+⋯,asϱ⟶0.We, therefore, are able to find the solution of equation (46) by solving leading order equation (47) subject to the inner boundary condition(63)ψ0′=Ks∗ψ0,onϱ⟶0.It is now expected that equation (47) subject to the inner boundary conditions (63) has multiple solutions. The hysteresis bifurcated point can be determined in a similar way. We obtain the critical values Ks∗H=0.129 and KH=4.712 for ϕ=0.0, Ks∗H=0.102 and KH=4.311 for ϕ=0.1, and Ks∗H=0.080 and KH=3.881 for ϕ=0.2, which further give(64a)KH⟶4.712,KsH∼0.129ScA1/2,asScA⟶0,whenϕ=0.0,(64b)KH⟶4.311,KsH∼0.102ScA1/2,asScA⟶0,whenϕ=0.1,(64c)KH⟶3.881,KsH∼0.080ScA1/2,asScA⟶0,whenϕ=0.2.


## 6. Conclusions

The homogeneous-heterogeneous reactions in a nanofluid with the suspended single-walled carbon nanotubes at the stagnation-point region of a plane surface have been investigated. The homogeneous (bulk) reaction is assumed to be isothermal cubic autocatalytic, while the heterogeneous (surface) reaction is single, isothermal, and first order. The homogeneous nanofluid model is employed with equal diffusion coefficients being considered for both reactants and autocatalysts. By doing similarity transformation, multiple solutions of the model have been captured for particular sets of physical parameters, which are overlooked in previous publications on modelling homogeneous-heterogeneous reactions. Furthermore, we have discovered the limitation of Magyari’s conclusion [25] about that the solutions of the homogeneous nanofluid model can be recovered from those of the regular fluid model. This statement is only valid for a few special situations.In summary, the main novel aspects of this study are as follows:(i)
Multiple solutions for modelling homogeneous-heterogeneous reactions are found for some sets of physical parameters, and their origins are discussed(ii)
Homogeneous and heterogeneous reaction ratesK and Ks, which are the obvious key factors to affect the chemical reactions, are important parameters to give multiple solutions(iii)
Nanoparticle volume fractionϕ plays a significant role on altering solution stability, and it also can affect the stability of chemical reactions(iv)
Magyari’s conclusion [25] on that the solutions of the homogeneous nanofluid model can be recovered from those of the regular fluid model is found to be limited and only valid for a few special cases



---

*Source: 3265143-2020-04-30.xml*