A Mathematical Study for Three-Dimensional Boundary Layer Flow of Jeffrey Nanofluid

1 Introduction

The suspension of ultrafine nanoparticles (<100 nm) in a base fluid is known as nanofluid. Nanofluid is a new variety of heat transfer fluids. To achieve the cooling rate requirements in industry, the thermal performance of ordinary heat transfer fluids are not suitable because such fluids have lower thermal conductivity. The concept of insertion of nanometer-sized metallic particles in fluid leads to an increase in the thermal conductivity of base liquids. In fact, the thermal conductivity of the base fluid is dramatically enhanced due to the presence of nanoparticles. The nanofluids are commonly used to explore the properties of colloidal suspension. Such fluids are a very interesting alternative in micro-electromechanical systems, nuclear reactors, electronic cooling equipment, transportation, heating and cooling process of energy conversion, etc. At present, the analysis of nanofluids is a topic of great interest to modern researchers due to its importance in industry and medicine. In addition, the magneto nanofluids have great importance in the construction of power generators, petroleum reservoirs, cancer therapy, wound treatment, gastric medications, elimination of tumours with hyperthermia, asthma treatment, and sterilised devices.

The variations in the thermal conductivities and viscosities of liquids through the dispersion of ultra-fine particles in the base fluids have been studied by Masuda et al. [1]. Choi [2] showed that the insertion of nanoparticles in base fluid is the best candidate to enhance the thermal performance of base fluids. Then Buongiorno [3] developed a mathematical model to explore the thermal properties of base fluids. This model involves the Brownian motion and thermophoresis effects. Afterwards the researchers investigated the flows of nanofluids under different conditions and for types of nanoparticles. Khan and Pop [4] used the Buongiorno model [3] to investigate the boundary layer flow of nanofluid past a stretching sheet. The Keller box method is used for the numerical solutions of the modelled differential system. Boundary layer flow of viscous nanofluid past a stretching surface with convective boundary condition was addressed by Makinde and Aziz [5]. They showed that the convective boundary condition is more realistic than the constant surface temperature or heat flux conditions. Mustafa et al. [6] explored the characteristics of thermophoresis and Brownian motion in stagnation point flow of viscous fluid past a stretching surface. Rashidi et al. [7] carried out a study to simulate the axisymmetric mixed convection nanofluid flow by a vertical cylinder. Further investigations on nanofluids can be seen in further studies [8–20]. Very recently, Kuznetsov and Nield [21] provided the revised model of natural convective boundary-layer flow of nanofluid past a vertical plate subject to the new proposed boundary condition. The authors claimed that the newly proposed boundary condition in study [21] is more realistic than the previous models. This condition describes that the nanofluid particle fraction on the boundary is passively rather than actively controlled, i.e., it is no longer assumed that one can control the value of the nanoparticle fraction at the wall but rather that the nanoparticle flux at the wall is zero. This change necessitates a rescaling of the parameters that are involved [21].

Boundary layer flow of non-Newtonian fluids over a stretching surface has received special attention from researchers. This is because of its occurrence in the industrial and engineering applications such as in hot rolling, wire drawing, manufacturing of plastic and rubber sheets, annealing and thinning of copper wires, continuous cooling of fiber spinning, drawing of stretching sheets through quiescent fluids, the boundary layer along a liquid film condensation process, and aerodynamic extrusion of plastic films. Several fluids, such as paints, paper pulp, shampoos, ketchup, apple sauce, certain oils, and polymer solutions are examples of non-Newtonian fluids. The non-Newtonian fluids are further divided into three categories namely; differential, rate, and integral types. The Jeffrey fluid model [22–25] is the subclass of rate type fluids. This fluid model describes the properties of ratio of relaxation to retardation times and retardation time.

The objective of the present study is to develop a mathematical model for magnetohydrodynamic (MHD) three-dimensional boundary layer flow of Jeffrey nanofluid over a stretching surface. Effects of thermophoresis and Brownian motion are introduced. To our knowledge, no such analysis of Jeffrey fluid has been performed in the literature yet. Therefore, here we prefer to utilise the recent condition studied by Kuznetsov and Nield [21] in the three-dimensional flow of Jeffrey fluid. Similarity variables are employed to convert the nonlinear partial differential system into an ordinary differential system. The resulting strong nonlinear system is computed and analysed through the series solutions by the homotopy analysis method (HAM) [26–30]. Convergence analysis of the obtained series solutions is made graphically and numerically. The discussion of plots and numerical values with respect to various parameters of interest is examined.

2 Mathematical Modelling

Let us consider the steady three-dimensional flow of an incompressible Jeffrey nanofluid. The flow is caused by a bidirectional stretching surface. The fluid is considered electrically conducting in the presence of constant magnetic field B₀ applied in the z-direction. In addition, the Hall and electric field effects are ignored. Induced magnetic field is not considered for small magnetic Reynolds numbers. Brownian motion and thermophoresis effects are taken into account. We adopt the Cartesian coordinate system in such a way that x- and y-axes are taken along the stretching surface and z-axis is normal to it. The sheet at z= 0 is stretched in the x- and y-directions with velocities U_w and V_w, respectively. The governing boundary layer expressions for the Jeffrey fluid are written as follows:

(1)∂u∂x + ∂v∂y + ∂w∂z = 0, (1)

(2)u∂u∂x + v∂u∂y + w∂u∂z  = ν1 + λ1(∂2u∂z2 + λ2(u∂3u∂x∂z2 + v∂3u∂y∂z2 + w∂3u∂z3+∂u∂z∂2u∂x∂z + ∂v∂z∂2u∂y∂z +∂w∂z∂2u∂z2 )) − σB02ρfu, (2)

(3)u∂v∂x + v∂v∂y + w∂v∂z  = ν1 + λ1(∂2v∂z2 + λ2(u∂3v∂x∂z2 + v∂3v∂y∂z2 + w∂3v∂z3+∂u∂z∂2v∂x∂z + ∂v∂z∂2v∂y∂z + ∂w∂z∂2v∂z2)) − σB02ρfv, (3)

(4)u∂T∂x + v∂T∂y + w∂T∂z = α∂2T∂z2 + (ρc)p(ρc)f(DB(∂T∂z∂C∂z) + DTT∞(∂T∂z)2), (4)

(5)u∂C∂x + v∂C∂y + w∂C∂z = DB(∂2C∂z2) + DTT∞(∂2T∂z2). (5)

The boundary conditions for the present flow analysis are [21, 23]

(6)u = Uw(x) = ax, v = Vw(y) = by, w = 0, T = Tw(x), DB∂C∂z + DTT∞∂T∂z = 0 at z = 0, (6)

(7)u → 0, v → 0, T → T∞, C → C∞   as z → ∞, (7)

where u, v, and w are the velocity components in the x-, y-, and z-directions, respectively, v(=μ/ρ_f) the kinematic viscosity, μ the dynamic viscosity, ρ_f the density of base fluid, λ₁ the ratio of relaxation and retardation times, λ₂ the retardation time, σ the electrical conductivity, T the temperature, α=k/(ρc)_f the thermal diffusivity of fluid, k the thermal conductivity, (ρc)_f the heat capacity of the fluid, (ρc)_p the effective heat capacity of nanoparticles, D_B the Brownian diffusion coefficient, C the nanoparticle concentration, D_T the thermophoretic diffusion coefficient, T_w and T_∞ the temperatures of the surface and far away from the surface, and C_∞ the nanoparticle concentration far away from the surface. The subscript w denotes the wall condition. Here we assume that the surface stretching velocities and wall temperature are

(8)Uw(x) = ax, Vw(y) = by, Tw(x) = T∞ + T0x. (8)

where a, b and T₀ are the positive constants.

The dimensionless variables can be defined as follows:

(9)u = axf′(η),  v = ayg′(η), w = −(aν)1/2(f(η) + g(η)),θ(η) = T − T∞Tw − T∞,  ϕ(η) = CC∞, η = (aν)1/2z. (9)

Equation (1) is automatically satisfied and (2)–(8) have the following forms:

(10)f‴ + (1 + λ1) ((f + g)f″ − f′2) + β(f″2 − (f + g)fiv − g′f‴)  − (1 + λ1)M2f′ = 0, (10)

(11)g‴ + (1 + λ1) ((f + g)g″ − g′2) + β(g″2 − (f + g) giv − f′g‴) − (1 + λ1) M2g′ = 0, (11)

(12)θ″ + Pr((f + g)θ′ − f′θ + Nbθ′ϕ′ + Ntθ′2) = 0, (12)

(13)ϕ″ + LePr(f + g)ϕ′ + NtNbθ″ = 0, (13)

(14)f = 0, g = 0, f′ = 1, g′ = c, θ = 1,  Nbϕ′+Ntθ′ = 0 at η = 0, (14)

(15)f′ → 0, g′ → 0, θ → 0, ϕ → 0 as η → ∞. (15)

In the previous expressions β is the Deborah number, M is the magnetic parameter, c is the ratio of stretching rates, Pr is the Prandtl number, Nb is the Brownian motion parameter, Nt is the thermophoresis parameter, Le is the Lewis number, and prime stands for differentiation with respect to η. These parameters can be expressed by the following definitions:

(16)β = λ2a,  M2 = σB02ρfa,  c = ba,  Pr = να,Nb = (ρc)pDBC∞(ρc)fν,  Nt = (ρc)pDT(Tw − T∞)(ρc)fνT∞,  Le = αDB.} (16)

The local Nusselt number Nu_x is defined as

(17)Nux = −x(Tw − T∞)∂T∂z|z = 0 = −(Rex)1/2θ′(0). (17)

It is noted that the dimensionless mass flux represented by a Sherwood number Sh_x is now identically zero and Re_x=U_wx/v is the local Reynolds number.

3 Development of Series Solutions

The initial guesses and the linear operators for homotopic solutions are selected as follows:

(18)f0(η) = 1 − e−η, g0(η) = c(1 − e−η), θ0(η) = e−η, ϕ0(η) = −NtNbe−η, (18)

(19)Lf=f‴−f′, Lg= g‴−g′, Lθ=θ″−θ, Lϕ=ϕ″−ϕ. (19)

The above operators satisfy the properties given as follows:

(20)Lf[C1 + C2eη + C3e−η] = 0, Lg[C4 + C5eη + C6e−η] = 0,Lθ[C7eη + C8e−η] = 0, Lϕ[C9eη + C10e−η] = 0,} (20)

in which C_i (i= 1 – 10) elucidate the arbitrary constants.

We can define the following zero-order deformation problems [27, 28] as

(21)(1 − p)Lf[f^(η, p) − f0(η)] = pℏfNf[f^(η, p), g^(η, p)], (21)

(22)(1 − p)Lg[g^(η, p) − g0(η)] = pℏgNg[f^(η, p), g^(η, p)], (22)

(23)(1 − p)Lθ[θ^(η, p) − θ0(η)] = pℏθNθ[f^(η, p), g^(η, p), θ^(η, p), ϕ^(η, p)], (23)

(24)(1 − p)Lϕ[ϕ^(η, p) − ϕ0(η)] = pℏϕNϕ[f^(η, p), g^(η, p), θ^(η, p), ϕ^(η, p)], (24)

(25)f^(0, p) = 0,  f^′(0, p) = 1,  f^′(∞, p) = 0, g^(0, p) = 0, g^′(0, p) = c, g^′(∞, p) = 0,θ^(0, p) = 1,  θ^(∞, p) = 0, Nbϕ^′(0, p)+Ntθ^′(0, p) = 0, ϕ^(∞, p) = 0,} (25)

(26)Nf[f^(η; p), g^(η; p)] = ∂3f^∂η3 + (1 + λ1) ((f^ + g^) ∂2f^∂η2 − (∂f^∂η)2)+β((∂2f^∂η2)2 − (f^ + g^) ∂4f^∂η4 − ∂g^∂η∂3f^∂η3) − (1 + λ1)M2∂f^∂η, (26)

(27)Ng[f^(η; p), g^(η; p)] = ∂3g^∂η3 + (1 + λ1)((f^ + g^) ∂2g^∂η2 − (∂g^∂η)2)+β((∂2g^∂η2)2 − (f^ + g^) ∂4g^∂η4 − ∂f^∂η∂3g^∂η3)− (1 + λ1) M2∂g^∂η, (27)

(28)Nθ[f^(η; p), g^(η; p), θ^(η, p), ϕ^(η, p)]  = ∂2θ^∂η2 + Pr(f^ + g^) ∂θ^∂η − Pr(∂f^∂η)θ^+NbPr∂θ^∂η∂ϕ^∂η   + NtPr(∂θ^∂η)2, (28)

(29)Nϕ[f^(η; p), g^(η; p), θ^(η, p), ϕ^(η, p)]  = ∂2ϕ^∂η2 + LePr(f^ + g^) ∂ϕ^∂η + NtNb∂2θ^∂η2. (29)

In the previous expressions p denotes the embedding parameter, ℏ_f, ℏ_g, ℏ_θ and ℏ_ϕ the nonzero auxiliary parameters and N_f, N_g, N_θ and N_ϕ the nonlinear operators. Setting p= 0 and p= 1 we have

(30)f^(η; 0) = f0(η),  f^(η; 1) = f(η), (30)

(31)g^(η; 0) = g0(η),  g^(η; 1) = g(η), (31)

(32)θ^(η, 0) = θ0(η),  θ^(η, 1) = θ(η), (32)

(33)ϕ^(η, 0) = ϕ0(η),  ϕ^(η, 1) = ϕ(η). (33)

When p varies from 0 to 1 then f^(η; p),g^(η; p),θ^(η, p) and ϕ^(η, p) vary from the initial guesses f₀(η), g₀(η), θ₀(η) and ϕ₀(η) to the final solutions f(η), g(η), θ(η) and ϕ(η), respectively. Taylor series expansion gives

(34)f^(η; p) = f0(η) + ∑m = 1∞fm(η)pm,  fm(η) = 1m!∂mf^(η, p)∂pm|p = 0, (34)

(35)g^(η; p) = g0(η) + ∑m = 1∞gm(η)pm,  gm(η) = 1m!∂mg^(η, p)∂pm|p=0, (35)

(36)θ^(η, p) = θ0(η)+∑m = 1∞θm(η)pm,  θm(η) = 1m!∂mθ^(η, p)∂pm|p = 0, (36)

(37)ϕ^(η, p) = ϕ0(η) + ∑m = 1∞ϕm(η)pm, ϕm(η) = 1m!∂mϕ^(η, p)∂pm|p = 0. (37)

The convergence of (34)–(37) strongly depends upon the suitable choices of ℏ_f, ℏ_g, ℏ_θ, and ℏ_ϕ. Considering that ℏ_f, ℏ_g, ℏ_θ, and ℏ_ϕ are chosen in such a manner that (34)–(37) converge at p= 1 then [29, 30]:

(38)f^(η; p) = f0(η) + ∑m = 1∞fm(η), (38)

(39)g^(η; p) = g0(η) + ∑m = 1∞gm(η), (39)

(40)θ^(η, p) = θ0(η) + ∑m = 1∞θm(η), (40)

(41)ϕ^(η, p) = ϕ0(η) + ∑m = 1∞ϕm(η). (41)

The general expressions of solutions can be expressed as follows which involve fm∗(η),gm∗(η),θm∗(η) and ϕm∗(η) as the special functions.

(42)fm(η) = fm∗(η) + C1 + C2eη + C3e−η, (42)

(43)gm(η) = gm∗(η) + C4 + C5eη + C6e−η, (43)

(44)θm(η) = θm∗(η) + C7eη + C8e−η, (44)

(45)ϕm(η) = ϕm∗(η) + C9eη + C10e−η. (45)

4 Convergence Analysis

The series solutions (38)–(41) involve the auxiliary parameters ℏ_f,ℏ_g,ℏ_θ, and ℏ_ϕ. These parameters are useful in adjusting and controlling the convergence of the series solutions. The proper values of these parameters play a vital role in constructing the convergent solutions through homotopy analysis method (HAM). To choose the suitable values of ℏ_f, ℏ_g, ℏ_θ, and ℏ_ϕ, the ℏ– curves are drawn at 20^th order of approximations. Figures 1 and 2 clearly show that the convergence region lies within the domain –1.35 ≤ ℏ_f ≤ –0.15, –1.45 ≤ ℏ_g ≤ –0.10, –1.40 ≤ ℏ_θ ≤ –0.45, and –1.40 ≤ ℏ_ϕ ≤ –0.25. Further, the presented solutions are convergent in the whole domain of η when ℏ_f=ℏ_g= –1.0 =ℏ_θ=ℏ_ϕ. Table 1 shows that the 34^th order of approximations are sufficient for the convergent series solutions.

Figure 1:

The h – curves for the functions f (η) and g(η).

Figure 2:

The h – curves for the functions θ(η) and ϕ(η).

5 Discussion

The purpose of this section is to explore the contributions of various influential parameters including Deborah number β, ratio of relaxation to retardation times λ₁, magnetic parameter M, ratio parameter c, Lewis number Le, Prandtl number Pr, Brownian motion parameter Nb, and thermophoresis parameter Nt on the dimensionless temperature and nanoparticle concentration distributions. Figures 3–16 are plotted for such a purpose.

Figure 3:

Influence of λ₁ on θ(η).

Figure 4:

Influence of β on θ(η).

Figure 5:

Influence of M on θ(η).

Figure 6:

Influence of c on θ(η).

Figure 7:

Influence of Nt on θ(η).

Figure 8:

Influence of Pr on θ(η).

Figure 9:

Influence of λ₁ on ϕ(η).

Figure 10:

Influence of β on ϕ(η).

Figure 11:

Influence of M on ϕ(η).

Figure 12:

Influence of c on ϕ(η).

Figure 13:

Influence of Le on ϕ(η).

Figure 14:

Influence of Pr on ϕ(η).

Figure 15:

Influence of Nt on ϕ(η).

Figure 16:

Influence of Nb on ϕ(η).

Figure 3 shows that an increase in ratio of relaxation to retardation times creates an enhancement in the temperature θ(η) and thermal boundary layer thickness. The case λ₁= 0 corresponds to minimum temperature and thinner thermal boundary layer thickness. Physically, an increase in λ₁ leads to an increase in relaxation time and a decrease in retardation time due to which higher temperature and thicker thermal boundary layer thickness. The variation in temperature profile θ(η) due to the Deborah number β is shown in Figure 4. It is noticed from Figure 4 that an increase in β shows a reduction in the temperature profile θ(η) and the related thermal boundary layer thickness. It can be analysed from Figure 5 that the temperature θ(η) and thermal boundary layer thickness are higher for the larger values of the magnetic parameter. Here M>0 is for hydromagnetic flow and M= 0 corresponds to hydrodynamic flow situation. We observed that the temperature and thermal boundary layer thickness are higher for hydromagnetic flow in comparison to the hydrodynamic case. The magnetic parameter depends upon the Lorentz force and larger magnetic parameters have a stronger Lorentz force that caused a tendency to increase the temperature in nanofluid motion. Influence of ratio parameter on the temperature profile θ(η) is plotted in Figure 6. From this figure it is noticed that the temperature and thermal boundary layer thickness are decreasing functions of the ratio parameter. For c= 0, a two-dimensional flow situation is recovered. Here we observed that the thermal boundary layer thickness is more for a two-dimensional case in comparison to three-dimensional flow.

Thermophoresis parameter Nt is a key parameter for analysing the temperature profile in nanofluid flow. The influence of thermophoresis parameter Nt on the temperature profile θ(η) is presented in Figure 7. The temperature of fluid increases when we increase the values of Nt. An increase in Nt produces an enhancement in the thermophoresis force which tends to move nanoparticles from hot to cold areas and, consequently, it enhances the temperature θ(η) and related thermal boundary layer thickness. The impact of the Prandtl number Pr on the temperature θ(η) is depicted in Figure 8. Increment in the Prandtl number creates a major effect on the temperature θ(η) and thermal boundary layer thickness. The thermal boundary layer thickness reduces with a larger Prandtl number and it occurs due to the decrease of thermal diffusivity for the higher values of a Prandtl number.

Figure 9 clearly shows that nanoparticle concentration ϕ(η) and its related boundary layer thickness are increasing functions of ratio of relaxation to retardation times. A comparison of Figures 3 and 9 clearly depicts that the impacts of ratio of relaxation to retardation times on the temperature and nanoparticle concentration are very similar. Influence of Deborah number β on the nanoparticle concentration ϕ(η) is plotted in Figure 10. Here the nanoparticle concentration and its related boundary layer thickness are reduced with the increase in Deborah number. Nanoparticle concentration ϕ(η) is an increasing function of the magnetic parameter (see Fig. 11). The nanoparticle concentration ϕ(η) enhances when the magnetic parameter increases.

The change in nanoparticle concentration ϕ(η) corresponding to different values of ratio parameter is examined in Figure 12. We have noticed that larger values of ratio parameter show a reduction in the nanoparticle concentration and its related boundary layer thickness. Figure 13 shows that nanoparticle concentration is reduced with an increase in the values of the Lewis number. The Lewis number depends upon the Brownian diffusion coefficient. Higher Lewis numbers lead to lower Brownian diffusion coefficients which shows a weaker nanoparticle concentration. Influence of the Prandtl number Pr on the nanoparticle concentration is depicted in Figure 14. Increase in the Prandtl number results in a reduction of the nanoparticle concentration. The nanoparticle concentration exhibits overshoot near the surface for larger values of Pr, though the nanoparticle volume boundary layer thickness reduces. The effect of thermophoresis parameter Nt on the nanoparticle concentration ϕ(η) is displayed in Figure 15. The nanoparticle concentration and its related boundary layer thickness are increasing functions of Nt. Increase in Nt creates an enhancement in the thermophoresis force which tends to increase the nanoparticle concentration and its related boundary layer thickness. The impact of Brownian motion parameter Nb on the nanoparticle concentration ϕ(η) is plotted in Figure 16. The nanoparticle concentration decreases with the increasing values of Nb. Through nanoparticles, the Brownian motion in the system takes place and with the increase in Nb, the Brownian motion is affected, and consequently the nanoparticle concentration boundary layer thickness decreases.

Table 2 is computed to investigate the heat transfer rate at the wall (local Nusselt number) for different values of β, λ_1,M, c, Nt,Nb, Le, and Pr. The heat transfer rate at the wall is constant for the variations in Brownian motion parameter. Effects of Lewis and Prandtl numbers on heat transfer rate are very opposite.

6 Concluding Remarks

Magnetohydrodynamic (MHD) three-dimensional boundary layer flow of Jeffrey nanofluid by a bidirectional stretching surface is addressed. The main points of present analysis are listed following:

• Temperature and nanoparticle concentration are enhanced with an increase in the ratio of relaxation to retardation times λ_1.

• Larger values of Deborah number β leads to a reduction in the temperature and nanoparticle concentration.

• Temperature θ(η) and thermal boundary layer thickness are decreasing functions of Prandtl number Pr.

• Temperature field increases when the values of Nt are increased.

• An increase in Lewis and Prandtl numbers show a decrease in nanoparticle concentration ϕ(η).

• Nanoparticle concentration ϕ(η) is reduced for an increase in Nb while it is increased with larger Nt.

• The heat transfer rate at the wall is constant for Brownian motion parameter Nb but it is a decreasing function of Nt.