1 Introduction

The complex nature of non-Newtonian fluids has posed interesting mathematical challenges for mathematicians, engineers, and researchers. This is because that non-Newtonian fluids play a vital role in the applications of physiology, biology, and industry. Common examples of such a type of applications include blood, shampoos, sauces, drilling muds, certain oils, lubricants, polymer solutions, and colloidal suspensions. The simple constitutive expression of Navier-Stokes is incapable to predict the mechanism of all such materials. This fact leads to the development of various non-Newtonian fluid models according to their physical nature. Here, we adopt the Oldroyd-B fluid model, which falls into the category of the rate type non-Newtonian liquids. The main feature of the Oldroyd-B fluid is to characterize the nature of the stress relaxation and retardation, which cannot be explored by the Maxwell fluid model. The idea of the boundary layer flow of the Oldroyd-B liquid induced by the linear stretching of a surface was first initiated by Sajid et al.^[1]. They reported the two-dimensional boundary layer flow of an Oldroyd-B liquid near a stagnation point, and developed the numerical solutions of the governing flow expressions. Many investigators have extended the work of Sajid et al.^[1]. Shehzad et al.^[2] discussed the effects of the temperature-dependent thermal conductivity on the three-dimensional flow of the Oldroyd-B fluid model, and considered that the flow generation was due to the bidirectional stretching of the surface. Hayat et al.^[3] addressed the effects of the temperature stratification in a steady-state stagnation point of the Oldroyd-B liquid with mixed convection, and developed the homotopic expressions of the solutions for velocity and temperature. Motsa et al.^[4] developed a spectral relaxation technique for the temperature-dependent three-dimensional flow of the Oldroyd-B liquid with heat source/sink effects. Abbasi et al.^[5] discussed the Cattaneo-Christov heat flux theory for the steady flow of the Oldroyd-B fluid over a moving sheet, and obtained the velocity and temperature by employing the optimal homotopic algorithm.

The thermal convection problems of fluid flows are very prominent in a number of industrial, engineering, and energy storage processes. Sheikholeslami et al.^[6] employed the lattice-Boltzmann technique to analyze the problem of the natural convection flow of a viscous nanoliquid under applied magnetic field. They considered the Cu-water nanoparticles filled in an annulus. Sheikholeslami et al.^[7] analyzed the forced convection non-uniform magnetohydrodynamic flow of a nanoliquid in a lid driven annulus. Mahanthesh et al.^[8] studied the mixed convective squeezing flow of a three-dimensional viscous nanofluid filled in a rotating channel, and presented numerical computations to examine the effects of various pertinent parameters. Rashidi et al.^[9] discussed the effects of mixed convection on the hydromagnetic flow of the Al₂O₃-water nanoliquid induced by a channel with sinusoidal walls and heat transfer. Abbasi et al.^[10] developed the homotopic algorithm to analyze the effects of double stratification in the mixed convective flow of the Maxwell nanoliquid over a moving surface with heat generation/absorption. Some recent investigations on convective flows can be found in Refs.[11]-[15].

Thermally radiative flows are generally encountered when the difference between the ambient temperature and the surface of the sheet is high. In several industrial processes, the thermal boundary layer thickness can be changed by use of thermal radiation. Examples of such industrial processes include nuclear reactors, power plants, satellites, missiles technology, gas turbines, etc. Abundant studies have been carried out in the literatures to predict the effects of thermal radiation (see Refs.[16]-[25] and the references therein). In these studies, the authors utilized the Rosseland approximations to linearize the thermal radiation term. In recent years, the investigation on non-linear thermal radiation has become a hot spot research topic. Cortell^[26] addressed the effects of the non-linear thermally radiative heat transfer in the steady laminar flow of a viscous liquid over a linear sheet. Mushtaq et al.^[27] analyzed the solar radiation effects in the two-dimensional flow of a viscous fluid, and presented a numerical analysis by taking the Brownian motion and thermophoretic effects into consideration. Shehzad et al.^[28] reported the non-linear thermal radiation effect in the three-dimensional Jeffrey nanoliquid over a bidirectional stretching surface. Hayat et al.^[29] reported the hydromagnetic three-dimensional viscoelastic fluid flow with non-linear thermal radiation. Mahanthesh et al.^[30] addressed the water-based nanofluid flow induced by the non-linear stretching surface under the effects of applied magnetic field and thermal radiation.

In this attempt, our main concern is to introduce the non-linear convection in the three-dimensional flow of an Oldroyd-B fluid induced by the stretching of the bidirectional stretching surface. We also consider the nonlinear thermal radiation and heat generation/absorption effects in the heat transfer expressions. The convective condition is employed at the boundary surface instead of the constant surface temperature condition. Different problems of the fluid flows have been treated by various numerical techniques^[31-40]. The present mathematical model is tackled through the fourth-and fifth-order formulae of the Runge-Kutta-Fehlberg techniques via the shooting algorithm. The results are plotted for multiple values of the dimensionless parameters to examine the curves of the velocities and temperature. The results are also discussed for the case of linear thermal radiation.

2 Flow and heat transfer analysis

The non-linear convection in an Oldroyd-B fluid past a stretching sheet is considered. A steady boundary layer flow is induced by the stretched surface at z=0, and it occupies the region z>0. The sheet is stretched in two directions with the velocities u_w=ax and v_w=by along the x-and y-directions correspondingly. We assume that T_f is the temperature of the convective surface, and T_∞ is the ambient fluid temperature. The magnetic Reynolds number is assumed to be so small that the induced magnetic field and the Hall current are negligible. The mathematical expressions of the conservation laws of mass, momentum, and energy subjected to the boundary layer assumptions are

(1)

(2)

(3)

(4)

where u, v, and w are, respectively, the velocity components along the x-, y-, and z-directions, T is the temperature, λ₁ is the relaxation time, λ₂ is the retardation time, ν is the kinematic viscosity, g is the acceleration due to gravity, β₀ and β₁ are the volumetric thermal expansion coefficients, α_m = k/(ρC_p) is the thermal diffusivity of the fluid, k is the thermal conductivity, ρ is the fluid density, C_p is the fluid specific heat, Q₀ is the heat generation/absorption coefficient, and q_r is the radiative heat flux. The present flow analysis is reduced to the Maxwell model by setting λ₂=0. Further, this analysis can be recovered for viscous liquids when λ₁=0=λ₂.

The thermal radiation heat flux expression through the Rosseland approximation is^[29]

(5)

where k₁ is the mean absorption coefficient, e_b = σT⁴ is the rate of the radiation emitted per square meter of surface, and σ is the Stefan-Boltzmann constant. The definition of e_b is through the Stefan-Boltzmann law, which states that all the objects with the temperature above absolute zero emit radiations at the rate proportional to the fourth power of its absolute temperature.

In view of Eq.(5), we have

(6)

The relevant boundary conditions for the present problem are

(7)

The governing partial differential equations suggest transformation into the corresponding nonlinear ordinary differential equations by choosing the following similarity variables^[28-29]:

(8)

where a prime denotes differentiation with respect to η. In view of the above relations, we obtain the following set of non-linear ordinary differential equations:

(9)

(10)

(11)

The boundary conditions for the present flow problem are

(12)

(13)

where the dimensionless parameters are

In the above equations, β₁ and β₂ are the Deborah numbers, c is the stretching ratio parameter, λ is the mixed convection parameter, Re_x is the Reynolds number, Gr_x is the Grashof number, γ is the non-linear convection parameter, R is the thermal radiation parameter, θ_w is the temperature ratio parameter, S is the heat source/sink parameter, Pr is the Prandtl number, and Bi is the Biot number.

The engineering interested physical quantity of the boundary value problems is the local Nusselt number Nu, which is defined by

(14)

where q_w is the surface heat flux. With the similarity variables, we obtain

(15)

Since the exact solution of the complicated nonlinear boundary value problem presented by Eqs.(9)-(13) is impracticable, we intend to handle this problem numerically.

3 Numerical analysis

The nonlinear boundary value problem is solved numerically by use of the Runge-Kutta-Fehlberg method along with the shooting technique. First, the non-linear differential equations are discretized to ten first-order linear differential equations. Then, the unknown initial conditions are calculated by use of the iterative technique, i.e., the shooting method, with some appropriate initial guesses. The fourth-and fifth-order formulae of the Runge-Kutta-Fehlberg method are^{[18, 30]}

(16)

(17)

where

The inner iteration is counted until the nonlinear solution converges with a convergence criterion of 10^-6. In addition, the step size is set to be ∆η=0.001. Another challenge to solve the system is fixing the appropriate finite values of η_∞. In this study, the asymptotic boundary conditions are replaced by η₈ in such a way that f ′(8) = g′(8) = f "(8) = g"(8) = θ(8) = 0. This ensures that all numerical solutions approach the asymptotic values correctly. In order to check the accuracy of our numerical method, the values of -f"(0) and -g"(0) for different values of the stretching ratio parameter are compared with those of Ref.[41], where the numerical results are obtained by the homotopy perturbation method (HPM) and the obtained exact solutions are for Newtonian fluids, when λ=R=0. The results are presented in Table 1. From the table, we can see that the present solutions are in good agreement with those of Ref.[41] as a limiting case.

4 Discussion

The graphs of the velocity distributions f'(η), g'(η) and the temperature field θ(η) for multiple values of the Deborah numbers β₁ and β₂, the ratio parameter c, the mixed convective parameter λ, the non-linear convection parameter γ, the radiation parameter R, the heat source/sink parameter S, the Biot number Bi, and the temperature ratio θ_w are visualized in Figs. 1-9.

Figure 1 is presented to explore the effects of the Deborah number β₁ on the velocities f'(η), g'(η) and the temperature θ (η). From the figure, we can see that, when β₁ increases, the velocities decrease first, then increase when β₁ is large enough. This is because that, the Deborah number β₁ depends on the relaxation time, and the relaxation time enhances for higher β₁. Such an enhancement in the relaxation time leads to lower velocities and higher temperature. It is also observed that the values of f'(η) at the wall are higher than those of g'(η) and θ(η) at the wall. The effects of β₂ on f'(η), g'(η), and θ (η) are visualized in Fig. 2.

From this figure, we see that f'(η) and g'(η) increase when β₂ increases. The variations in g'(η) are quite prominent in comparison with the changes in the curves of f'(η). Moreover, an increase in β₂ leads to lower temperature and thinner thermal boundary layer. From the definition of β₂, we can see that the retardation factor is higher for larger β₂, which may create decreases in the temperature and the thermal boundary layer. The present results can be modified to the results of the Maxwell fluid by taking β₂ =0. The analysis of the three-dimensional flow of the viscous liquid can be retrieved by setting β₁ =0=β₂.

The velocity f'(η) and the temperature θ (η) decay while the velocity g'(η) increases remarkably when c increases (see Fig. 3). It is due to the fact that an increase in c from zero leads to the movement of the lateral surface in the y-direction that corresponds to a higher velocity g'(η) and its associated boundary layer thickness. The present three-dimensional problem can be converted into a two-dimensional flow model when c=0. From Fig. 4, we can see that the velocity f'(η) increases remarkably while the velocity g'(η) and the temperature θ (η) decrease when the mixed convective parameter λ increases. This occurs due to the buoyancy force in λ. The curves of f'(η), g'(η), and θ (η) for various values of the non-linear convection parameter are given in Fig. 5. Here, the velocity f'(η) is an increasing function of the non-linear convection parameter. When γ increases, g'(η) and θ (η) decrease.

Figure 6 clearly shows that, when the radiative parameter increases, the values of f'(η) and θ (η) become higher, while the values of the velocity g'(η) become smaller. More heat is generated in the fluid due to the radiation that corresponds to the thicker momentum and the thermal boundary layer thicknesses. Similarly, an increase in the heat source/sink parameter leads to increases in the velocity f'(η) and the temperature θ (η) and a decrease in the velocity g'(η) (see Fig. 7). Here, S=0 implies no heat source/sink, and S>0 corresponds to heat source. The heat sink case occurs when the values of S are negative.

The variations in f'(η), g'(η), and θ (η) corresponding to different values of the Biot number Bi are explored in Fig. 8. The variations in the curves of g'(η) for multiple values of Bi are very small. The velocity f'(η) and the temperature θ (η) are enhanced significantly due to an increase in Bi. The heat transfer coefficient becomes stronger for larger Bi, which gives rise to the fluid velocity f'(η) and the temperature θ (η). The temperature ratio θ_w leads to remarkable changes in the curves of f'(η) and θ (η), while the profile of g'(η) changes very slowly (see Fig. 9).

The comparison of the present numerical results with those in Ref.[41] is made in Table 1. It is indicated that the present values of f"(0) and g"(0) for multiple values of c have an excellent match with the results of Ariel^[41]. Table 2 is presented to examine the values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0) for multiple values of β₁ and β₂ when c=0.6, Pr =1.2, S=0.3, R=0.4, θ_w =1.6, and Bi=0.4. Here, we have computed the values by considering λ =0 and γ =0.5, λ =γ=0.5, and λ =0.5 and γ=0. From Table 2, we can see that the values of f"(0) and -(1+Rθ_w³)θ '(0) when λ =0 and γ=0.5 or λ =0.5 and γ=0 are smaller than those when λ =γ=0.5, while the values of g"(0) when λ =0 and γ=0.5 or λ =0.5 and γ=0 are bigger than those when λ =γ=0.5. It can be also seen that the values of f"(0) and g"(0) decay with an enhancement in β₁ while boost up with an increase in β₂. The values of -(1+Rθ_w³)θ'(0) are enhanced for larger β₁, while decrease when β₂ increases.

The values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0) for different values of R by setting

are investigated in Table 3. In this Table, we make an analysis of the values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0) in absence of radiation, linear radiation, and non-linear radiation. The results show that the radiation term has no effect on the values of f"(0) and g"(0) when λ =0, and γ=0.5. The numerical values of -(1+Rθ_w³)θ'(0) are larger in the non-linear radiation situation than those in the cases of linear radiation and absence of radiation. In Table 4, we have studied the values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0) for various values of Bi by setting

Here, we notice that smaller values of the Biot number have greater effects on the values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0).

5 Conclusions

The role of non-linear convection and thermal radiation in the three-dimensional flow of the Oldroyd-B liquid is explored. The heat transfer phenomenon is examined under the heat source/sink and convective surface condition. Numerical computations have been carried out to analyze the solutions of the velocities and temperature. The results show that the Deborah numbers β₁ and β₂ have reverse effects on the velocities and temperature. It is also noted that the values of the velocity f'(η) at the wall are higher than the values of the velocity g'(η) and the temperature θ (η). Larger c leads to smaller f'(η) while bigger g'(η). The velocity f'(η) is an increasing function of the mixed convection parameter. When the non-linear convection parameter γ increases, the velocity g'(η) and the temperature θ (η) decrease. The temperature θ (η) increases significantly when Bi increases. The values of f"(0) and -(1+Rθ_w³)θ'(0) when λ=0 and γ=0.5 or λ=0.5 and γ=0 are smaller than those when λ=γ=0.5. The values of g"(0) when λ=0 and γ=0.5 or λ=0.5 and γ=0 are bigger than those when λ=γ=0.5. The results also show that the radiation term has no effect on the values of f"(0) and g"(0) when λ =0, and γ=0.5. The numerical values of -(1+Rθ_w³)θ'(0) in the case of non-linear radiation are larger than those in the cases of linear radiation and absence of radiation. It is noticed that smaller values of the Biot number have greater effects on the values of f"(0), g"(0), and -(1+Rθ_w³)θ'(0).