Unsteady MHD free convection boundary-layer flow of a nanofluid along a stretching sheet with thermal radiation and viscous dissipation effects

Abstract

In this work, we study the unsteady free convection boundary-layer flow of a nanofluid along a stretching sheet with thermal radiation in the presence of magnetic field. To obtain non-similar equations, continuity, momentum, energy, and concentration equations have been non-dimensionalized by usual transformation. The non-similar solutions are considered here which depend on the magnetic parameter M, radiation parameter R, Prandtl number P_r, Eckert number E_c, Lewis number L_e, Brownian motion parameter N_b, thermophoresis parameter N_t, and Grashof number G_r. The obtained equations have been solved by an explicit finite difference method with stability and convergence analysis. The velocity, temperature, and concentration profiles are discussed for different time steps and for the different values of the parameters of physical and engineering interest.

Background

Magnetohydrodynamics (MHD) boundary-layer flow of nanofluid and heat transfer over a linearly stretched surface have received a lot of attention in the field of several industrial, scientific, and engineering applications in recent years. Nanofluids have many applications in the industries since materials of nanometer size have unique chemical and physical properties. With regard to the sundry applications of nanofluids, the cooling applications of nanofluids include silicon mirror cooling, electronics cooling, vehicle cooling, transformer cooling, etc. This study is more important in industries such as hot rolling, melt spinning, extrusion, glass fiber production, wire drawing, and manufacture of plastic and rubber sheets, polymer sheet and filaments, etc.

Sakiadis [1] was the first author to analyze the boundary-layer flow on a continuous surface. Crane [2] obtained an exact solution of the boundary-layer flow of the Newtonian fluid caused by the stretching of an elastic sheet moving in its own plane linearly. Gorla et al. [3, 4] solved the non-similar problem of free convective heat transfer from a vertical plate embedded in a saturated porous medium with an arbitrarily varying surface temperature. Cheng and Minkowycz [5] also studied free convection from a vertical flat plate with applications to heat transfer from a dick.

Dissipation is the process of converting mechanical energy of downward-flowing water into thermal and acoustical energy. Various devices are designed in streambeds to reduce the kinetic energy of flowing waters, reducing their erosive potential on banks and river bottoms. Vajravelu and Hadjinicalaou [6] analyzed the heat transfer characteristics over a stretching surface with viscous dissipation in the presence of internal heat generation or absorption.

Takhar et al. [7] studied the radiation effects on the MHD free convection flow of a gas past a semi-infinite vertical plate. Ghaly [8] considered the thermal radiation effect on a steady flow, whereas Rapits and Massalas [9] and El-Aziz [10] analyzed the unsteady case. Sattar and Alam [11] presented unsteady free convection and mass transfer flow of a viscous, incompressible, and electrically conducting fluid past a moving infinite vertical porous plate with thermal diffusion effect. Na and Pop [12] analyzed an unsteady flow due to a stretching sheet. In the case of unsteady boundary-layer flow, Singh et al. [13] investigated the thermal radiation and magnetic field effects on an unsteady stretching permeable sheet in the presence of free stream velocity. The study of convective instability and heat transfer characteristics of nanofluids was considered by Kim et al. [14]. Jang and Choi [15] obtained nanofluid thermal conductivity and the effect of various parameters on it.

The natural convective boundary-layer flows of a nanofluid past a vertical plate have been described by Kuznetsov and Nield [16, 17]. In this model, the Brownian motion and thermophoresis are accounted with the simplest possible boundary conditions. They also studied the Cheng-Minkowycz problem for natural convective boundary-layer flow in a porous medium saturated by a nanofluid. Bachok et al. [18] have shown the steady boundary-layer flow of a nanofluid past a moving semi-infinite flat plate in a uniform free stream. It was assumed that the plate is moving in the same or opposite directions to the free stream to define the resulting system of non-linear ordinary differential equations.

Khan and Pop [19, 20] formulated the problem of laminar boundary-layer flow of a nanofluid past a stretching sheet. They also expressed free convection boundary-layer nanofluid flow past a horizontal flat plate. Hamad and Pop [21] discussed the boundary-layer flow near the stagnation-point flow on a permeable stretching sheet in a porous medium saturated with a nanofluid. Hamad et al. [22] investigated free convection flow of a nanofluid past a semi-infinite vertical flat plate with the influence of magnetic field. Very recently, Hady et al. [23] investigated the effects of thermal radiation on the viscous flow of a nanofluid and heat transfer over a non-linearly stretching sheet.

However, the aim of the present work is to study the unsteady free convection boundary-layer nanofluid flows along a stretching surface with the influence of magnetic field and radiation effect. An explicit finite difference procedure [24] has been taken to solve the obtained non-similar equations with stability and convergence analysis.

Methods

Presentation of the hypothesis

An unsteady two-dimensional MHD free convection laminar boundary-layer flow of a viscous incompressible and electrically conducting nanofluid along a vertical stretching sheet under the influence of thermal radiation and viscous dissipation is considered. The sketch of the physical configuration and coordinate system is shown in Figure 1. Introducing the Cartesian coordinate system, the x-axis is taken along the stretching sheet in the vertically upward direction, and the y-axis is taken as normal to the sheet. Two equal and opposite forces are introduced along the x-axis so that the sheet is stretched, keeping the origin fixed.

Figure 1

Physical model and coordinate system.

Instantaneously at time t > 0, the temperature of the plate and the species concentration are raised to T_w(> T _∞ ) and C_w(> C _∞ ), respectively, which are thereafter maintained constant, where T_w and C_w are the temperature and species concentration at the wall, respectively, and T_∞ and C_∞ are the temperature and species concentration far away from the plate, respectively.

A strong magnetic field is applied in the y direction. The uniform magnetic field strength (magnetic induction) B₀ can be taken as B = (0, B₀, 0). The Rosseland approximation is used to describe the radioactive heat flux q_r in the energy equation. Under the above assumptions and the usual boundary layer approximation, the MHD free convection unsteady nanofluid flow and heat and mass transfer with the radiation effect are governed by the following equations:

• The continuity equation

  $$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

  (1)

• The momentum equation

  $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\upsilon \frac{{\partial }^{2}u}{\partial y^2}+g\beta \left(T-T_{\infty }\right)-\frac{\sigma B_0^{2}u}{\rho }$$

  (2)

• The energy equation

  $$\begin{array}{cc}\frac{\partial T}{\partial t}\hfill & +u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{k}{\rho c_{\text{p}}}\frac{{\partial }^{2}T}{\partial y^2}\hfill \\ -\frac{1}{\rho c_{\text{p}}}\frac{\partial q_{\text{r}}}{\partial y}+\frac{\upsilon }{c_{\text{p}}}{\left(\frac{\partial u}{\partial y}\right)}^2\hfill \\ +\tau \left\{D_{\text{B}}\left(\frac{\partial T}{\partial y}·\frac{\partial C}{\partial y}\right)+\frac{D_{\text{T}}}{T_{\infty }}{\left(\frac{\partial T}{\partial y}\right)}^2\right\}\hfill \end{array}$$

  (3)

• The concentration equation

  $$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_{\text{B}}\frac{{\partial }^{2}C}{\partial y^2}+\frac{D_{\text{T}}}{T_{\infty }}\frac{{\partial }^{2}T}{\partial y^2}$$

  (4)

The initial and boundary conditions are

$$\begin{array}{l}t=0,u=Dx,v=0,T=T_{\infty },C=C_{\infty }\text{everywhere}\\ t\ge 0,u=0,v=0,T=T_{\infty },C=C_{\infty }\text{at}x=0\\ u=Dx,v=0,T=T_{\text{w}},C=C_{\text{w}}\text{at}y=0\\ u=0,v=0,T\to T_{\infty },C\to C_{\infty }\text{at}y\to \infty \text{,}\end{array}$$

(5)

where u and v are the velocity components in the x and y directions, respectively, v is the kinematic viscosity, k is the thermal conductivity, D_B is the Brownian diffusion coefficient, D_T is the thermophoresis diffusion coefficient, D(> 0) is the stretching constant, g is the acceleration due to gravity, ρ is the density of the fluid, and c_p is the specific heat at constant pressure.

The Rosseland approximation [25] is expressed for radiative heat flux and leads to the form

$$q_{\text{r}}=-\frac{4\sigma }{3{\kappa }^{*}}\frac{\partial T^4}{\partial y}\text{,}$$

(6)

where σ is the Stefan-Boltzmann constant and κ^* is the mean absorption coefficient. The temperature difference within the flow is sufficiently small such that T⁴ may be expressed as a linear function of the temperature, then Taylor's series for T⁴ is about T_∞ after neglecting higher order terms:

$$T^4=4{T_{\infty }}^3-3{T_{\infty }}^4\text{.}$$

(7)

Introducing the following non-dimensional variables, $X=\frac{x{U}_0}{\upsilon },Y=\frac{y{U}_0}{\upsilon },U=\frac{u}{U_0},V=\frac{v}{U_0}\text{,}$$\tau =\frac{t{U}_0^2}{\upsilon },\bar{T}=\frac{T-T_{\infty }}{T_{\text{w}}-T_{\infty }},\bar{C}=\frac{C-C_{\infty }}{C_{\text{w}}-C_{\infty }}\text{.}$

Then, Equations 1 to 5 become

$$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$$

(8)

$$\frac{\partial U}{\partial \tau }+U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=\frac{{\partial }^{2}U}{\partial Y^2}+G_{\text{r}}\bar{T}-MU$$

(9)

$$\begin{array}{cc}\frac{\partial \bar{T}}{\partial \tau }\hfill & +U\frac{\partial \bar{T}}{\partial X}+V\frac{\partial \bar{T}}{\partial Y}=\left[\left(\frac{1+R}{P_{\text{r}}}\right)·\left(\frac{{\partial }^2\bar{T}}{\partial Y^2}\right)\right]\hfill \\ +E_{\text{c}}{\left(\frac{\partial U}{\partial Y}\right)}^2+N_{\text{b}}\left(\frac{\partial \bar{T}}{\partial Y}·\frac{\partial \bar{C}}{\partial Y}\right)+N_{\text{t}}{\left(\frac{\partial \bar{T}}{\partial Y}\right)}^2\hfill \end{array}$$

(10)

$$\frac{\partial \bar{C}}{\partial \tau }+U\frac{\partial \bar{C}}{\partial X}+V\frac{\partial \bar{C}}{\partial Y}=\frac{1}{L_{\text{e}}}\left[\frac{{\partial }^2\bar{C}}{\partial Y^2}+\left(\frac{N_{\text{t}}}{N_{\text{b}}}\right)\frac{{\partial }^2\bar{T}}{\partial Y^2}\right]\text{.}$$

(11)

The non-dimensional boundary conditions are

$$\begin{array}{ll}\tau \le 0,U=0,V=0,\bar{T}=0,\bar{C}=0\hfill & \text{everywhere}\hfill \\ \tau >0,U=0,V=0,\bar{T}=0,\bar{C}=0\hfill & \text{at}X=0\hfill \end{array}$$

(12)

$$\begin{array}{ll}U=1,V=0,\bar{T}=1,\bar{C}=1\hfill & \text{at}Y=0\hfill \\ U=0,V=0,\bar{T}=0,\bar{C}=0\hfill & \text{as}Y\to \infty \text{,}\hfill \end{array}$$

(13)

where the magnetic parameter $M=\frac{\sigma B_0^2\nu }{\rho U_0^2}$, Grashof number $G_{\text{r}}=\frac{g\beta \left(T_{\text{w}}-T_{\infty }\right)\nu }{U_o^3}$, radiation parameter $R=\frac{16\sigma T_{\infty }^3}{3k{\kappa }^{*}}$, Prandtl number $P_{\text{r}}=\frac{\upsilon }{\alpha }$, Eckert number $E_{\text{c}}=\frac{U_0^2}{c_{\text{p}}\left(T_{\text{w}}-T_{\infty }\right)}$, Lewis number $L_{\text{e}}=\frac{\upsilon }{D_{\text{B}}}$, Brownian parameter $N_{\text{b}}=\frac{\tau D_{\text{B}}\left(C_{\text{w}}-C_{\infty }\right)}{\upsilon }$, and thermophoresis parameter $N_{\text{t}}=\frac{D_{\text{T}}}{T_{\infty }}\frac{\tau }{\nu }\left(T_{\text{w}}-T_{\infty }\right)$.

Methods: numerical computation

In order to solve the non-similar unsteady coupled non-linear partial differential equations (Equations 8, 9, 10, and 11), the explicit finite difference method has been developed. For this, a rectangular region of the flow field is chosen, and the region is divided into a grid of lines parallel to X and Y axes, where the X-axis is taken along the plate and the Y-axis is normal to the plate.

Here, the plate of height X_max(=100) is considered, i.e., X varies from 0 to 100 and assumed Y_max(=25) as corresponding to Y → ∞, i.e., Y varies from 0 to 25. There are m(=125) and n(=125) grid spacing in the X and Y directions, respectively, as shown in Figure 2. It is assumed that ΔX and ΔY are constant mesh sizes along the X and Y directions, respectively, and taken as follows: ΔX = 0.8(0 ≤ X ≤ 100) and ΔY = 0.2(0 ≤ Y ≤ 25) with the smaller timestep Δτ = 0.005.

Figure 2

Finite difference space grid.

Let U^′, V^′, ${\bar{T}}^{\prime }$, and ${\bar{C}}^{\prime }$ denote the values of U, V, $\bar{T}$, and $\bar{T}$ at the end of a time step, respectively. Using the explicit finite difference approximation, the following appropriate sets of finite difference equations are obtained:

$$\frac{{U^{\prime }}_{i,j}-{U^{\prime }}_{i-1,j}}{\Delta X}+\frac{V_{i,j}-V_{i,j-1}}{\Delta Y}=0$$

(14)

$$\begin{array}{c}\frac{{U^{\prime }}_{i,j}-U_{i,j}}{\Delta \tau }+U_{i,j}\frac{U_{i,j}-U_{i-1,j}}{\Delta X}+V_{i,j}\frac{U_{i,j+1}-U_{i,j}}{\Delta Y}\hfill \\ =\frac{U_{i,j+1}-2U_{i,j}+U_{i,j-1}}{{\left(\Delta Y\right)}^2}+G_{\text{r}}{T}^{\prime }-M{U}_{i,j}\hfill \end{array}$$

(15)

$$\begin{array}{c}\frac{{{\bar{T}}^{\prime }}_{i,j}-{\bar{T}}_{i,j}}{\Delta \tau }+U_{i,j}\frac{{\bar{T}}_{i,j}-{\bar{T}}_{i-1,j}}{\Delta X}+V_{i,j}\frac{{\bar{T}}_{i,j+1}-{\bar{T}}_{i,j}}{\Delta Y}\hfill \\ =\left(\frac{1+R}{P_{\text{r}}}\right)\left(\frac{{\bar{T}}_{i,j+1}-2{\bar{T}}_{i,j}+{\bar{T}}_{i,j-1}}{{\left(\Delta Y\right)}^2}\right)\hfill \\ +E_{\text{c}}{\left(\frac{U_{i,j+1}-U_{i,j}}{\Delta Y}\right)}^2\hfill \\ +N_{\text{b}}\left(\frac{{\bar{T}}_{i,j+1}-{\bar{T}}_{i,j}}{\Delta Y}·\frac{{\bar{C}}_{i,j+1}-{\bar{C}}_{i,j}}{\Delta Y}\right)\hfill \\ +N_{\text{t}}{\left(\frac{{\bar{T}}_{i,j+1}-{\bar{T}}_{i,j}}{\Delta Y}\right)}^2\hfill \end{array}$$

(16)

$$\begin{array}{c}\frac{{\overline{C^{\prime }}}_{i,j}-{\bar{C}}_{i,j}}{\Delta \tau }+U_{i,j}\frac{{\bar{C}}_{i,j}-{\bar{C}}_{i-1,j}}{\Delta X}+V_{i,j}\frac{{\bar{C}}_{i,j+1}-{\bar{C}}_{i,j}}{\Delta Y}\hfill \\ =\frac{1}{L_{\text{e}}}\left[\left(\frac{{\bar{C}}_{i,j+1}-2{\bar{C}}_{i,j}+{\bar{C}}_{i,j-1}}{{\left(\Delta Y\right)}^2}\right)\right.\hfill \\ \left.+\frac{N_{\text{t}}}{N_{\text{b}}}\left(\frac{{\bar{T}}_{i,j+1}-2{\bar{T}}_{i,j}+{\bar{T}}_{i,j-1}}{{\left(\Delta Y\right)}^2}\right)\right]\hfill \end{array}$$

(17)

with initial and boundary conditions

$$U_{i,j}^0=0,V_{i,j}^0=0,{\bar{T}}_{i,j}^0=0,{\bar{C}}_{i,j}^0=0$$

(18)

$$U_{0,j}^n=0,V_{0,j}^n=0,{\bar{T}}_{0,j}^0=0,{\bar{C}}_{0,j}^n=0$$

$$U_{i,0}^n=1,V_{i,0}^n=0,{\bar{T}}_{i,0}^n=1,{\bar{C}}_{i,0}^n=1$$

(19)

$$U_{i,L}^n=0,V_{i,L}^n=0,{\bar{T}}_{i,L}^n=0,{\bar{C}}_{i,L}^n=0,\text{where}L\to \infty \text{,}$$

where the subscripts i and j designate the grid points with X and Y coordinates, respectively, and the superscript n represents a value of time, τ = n · Δτ, where n = 0, 1, 2, ….

Stability and convergence analysis

Since an explicit procedure is being used, the analysis will remain incomplete unless the stability and convergence of the finite difference scheme is discussed. For the constant mesh sizes, the stability criteria of the scheme may be established as follows.

Equation 14 will be ignored since Δτ does not appear in it. The general terms of the Fourier expansion for U, $\bar{T}$, and $\bar{C}$ at a time arbitrarily called τ = 0 are all e^iαXe^iβY, apart from a constant, where $i=\sqrt{-1}$. At a time τ, these terms become

$$\begin{array}{c}\hfill U:\psi \left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\hfill \\ \hfill \bar{T}:\theta \left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\hfill \\ \hfill \bar{C}:\varphi \left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\text{,}\hfill \end{array}$$

(20)

and after the time step, these terms will become

$$\begin{array}{c}\hfill U:{\psi }^{\prime }\left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\hfill \\ \hfill \bar{T}:{\theta }^{\prime }\left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\hfill \\ \hfill \bar{C}:{\varphi }^{\prime }\left(\tau \right)e^{\mathit{i\alpha X}}{e}^{\mathit{i\beta Y}}\text{.}\hfill \end{array}$$

(21)

Substituting Equations 20 and 21 into Equations 15 to 17, with regard to the coefficients U and V as constants over any one time step, we obtain the following equations upon simplification:

$$\begin{array}{c}\frac{{\psi }^{\prime }\left(\tau \right)-\psi \left(\tau \right)}{\text{Δ}\tau }+U\frac{\psi \left(\tau \right)\left(1-e^{-i\alpha \text{Δ}X}\right)}{\text{Δ}X}\hfill \\ +V\frac{\psi \left(\tau \right)\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}=\frac{2\psi \left(\tau \right)\left(\text{cos}\beta \text{Δ}Y-1\right)}{{\left(\text{Δ}Y\right)}^2}\hfill \\ +G_{\text{r}}{\theta }^{\prime }\left(\tau \right)-M\psi \left(\tau \right)\hfill \end{array}$$

(22)

$$\begin{array}{c}\frac{{\theta }^{\prime }\left(\tau \right)-\theta \left(\tau \right)}{\text{Δ}\tau }+U\frac{\theta \left(\tau \right)\left(1-e^{-i\alpha \text{Δ}X}\right)}{\text{Δ}X}\hfill \\ +V\frac{\theta \left(\tau \right)\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}=\left(\frac{1+R}{P_{\text{r}}}\right)\frac{2\theta \left(\tau \right)\left(\text{cos}\beta \text{Δ}Y-1\right)}{{\left(\text{Δ}Y\right)}^2}\hfill \\ +E_{\text{c}}U\psi \left(\tau \right){\left\{\frac{\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}\right\}}^2\hfill \\ +N_{\text{b}}\bar{C}\theta \left(\tau \right){\left\{\frac{\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}\right\}}^2\hfill \\ +N_{\text{t}}\bar{T}\theta \left(\tau \right){\left\{\frac{\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}\right\}}^2\hfill \end{array}$$

(23)

$$\begin{array}{c}\frac{{\varphi }^{\prime }\left(\tau \right)-\varphi \left(\tau \right)}{\text{Δ}\tau }+U\frac{\varphi \left(\tau \right)\left(1-e^{-i\alpha \text{Δ}X}\right)}{\text{Δ}X}\hfill \\ +V\frac{\varphi \left(\tau \right)\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}=\frac{1}{L_{\text{e}}}\left[\left\{\frac{2\varphi \left(\tau \right)\left(\text{cos}\beta \text{Δ}Y-1\right)}{{\left(\text{Δ}Y\right)}^2}\right\}\right]\hfill \\ \left.+\left(\frac{N_{\text{t}}}{N_{\text{b}}}\right)·\left\{\frac{2\theta \left(\tau \right)\left(\text{cos}\beta \text{Δ}Y-1\right)}{{\left(\text{Δ}Y\right)}^2}\right\}\right]\hfill \end{array}$$

(24)

Equations 22, 23, and 24 can be written in the following forms:

$${\psi }^{\prime }=A\psi +G_{\text{r}}\text{Δ}\tau {\theta }^{/}$$

(25)

$${\theta }^{\prime }=B\theta +E\psi$$

(26)

$${\varphi }^{\prime }=J\varphi +K\theta \text{,}$$

(27)

where

$$\begin{array}{cc}A=\hfill & 1-U\frac{\Delta \tau }{\Delta X}\left(1-e^{-i\alpha \Delta X}\right)-V\frac{\Delta \tau }{\text{Δ}Y}\left(e^{i\beta \text{Δ}Y}-1\right)\hfill \\ +\frac{2\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\left(\text{cos}\beta \text{Δ}Y-1\right)-M\text{Δ}\tau \hfill \end{array}$$

$$\begin{array}{cc}B=\hfill & 1-U\frac{\text{Δ}\tau }{\text{Δ}X}\left(1-e^{-i\alpha \text{Δ}X}\right)-V\frac{\text{Δ}\tau }{\text{Δ}Y}\left(e^{i\beta \text{Δ}Y}-1\right)\hfill \\ +\left(\frac{1+R}{P_{\text{r}}}\right)\frac{2\left(\text{cos}\beta \text{Δ}Y-1\right)}{{\left(\text{Δ}Y\right)}^2}\text{Δ}\tau \hfill \\ +N_{\text{b}}\bar{C}{\left\{\frac{\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}\right\}}^2\text{Δ}\tau \hfill \\ +N_{\text{t}}\bar{T}{\left\{\frac{\left(e^{i\beta \text{Δ}Y}-1\right)}{\text{Δ}Y}\right\}}^2\text{Δ}\tau \text{,}\hfill \end{array}$$

$$E=E_{\text{c}}U\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}{\left(e^{i\beta \text{Δ}Y}-1\right)}^2\text{,}$$

$$\begin{array}{cc}J=\hfill & 1-U\frac{\text{Δ}\tau }{\text{Δ}X}\left(1-e^{-i\alpha \text{Δ}X}\right)-V\frac{\text{Δ}\tau }{\text{Δ}Y}\left(e^{i\beta \text{Δ}Y}-1\right)\hfill \\ +\frac{1}{L_{\text{e}}}\frac{2\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\left(\text{cos}\beta \text{Δ}Y-1\right)\text{,}\hfill \end{array}$$

and

$$K=\frac{1}{L_{\text{e}}}\left(\frac{N_{\text{t}}}{N_{\text{b}}}\right)\frac{2\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\left(\text{cos}\beta \text{Δ}Y-1\right)\text{.}$$

Again using Equation 26 in Equation 25,

$$\begin{array}{c}{\psi }^{\prime }=A\psi +G_{\text{r}}\text{Δ}\tau \left(B\theta +E\psi \right)\\ =C\psi +D\theta \text{,}\end{array}$$

where C = A + EG_rΔτ and D = BG_rΔτ.

Therefore, Equations 25, 26, and 27 can be expressed as

$${\psi }^{\prime }=C\psi +D\theta$$

(28)

$${\theta }^{\prime }=B\theta +E\psi$$

(29)

$${\varphi }^{\prime }=J\varphi +K\theta \text{.}$$

(30)

Hence, Equations 28, 29, and 30 can be expressed in a matrix notation, and these equations are

$$\left[\begin{array}{c}\hfill {\psi }^{\prime }\hfill \\ \hfill {\theta }^{\prime }\hfill \\ \hfill {\varphi }^{\prime }\hfill \end{array}\right]=\left[\begin{array}{c}\hfill CD0\hfill \\ \hfill 0B0\hfill \\ \hfill 0KJ\hfill \end{array},\right].\left[\begin{array}{c}\hfill \psi \hfill \\ \hfill \theta \hfill \\ \hfill \varphi \hfill \end{array}\right]\text{,}$$

(31)

that is, η^′ = Tη, where ${\eta }^{\prime }=\left[\begin{array}{c}\hfill {\psi }^{\prime }\hfill \\ \hfill {\theta }^{\prime }\hfill \\ \hfill {\varphi }^{\prime }\hfill \end{array}\right],T=\left[\begin{array}{c}\hfill CD0\hfill \\ \hfill 0B0\hfill \\ \hfill 0KJ\hfill \end{array},\right].\text{and}\eta =\left[\begin{array}{c}\hfill \psi \hfill \\ \hfill \theta \hfill \\ \hfill \varphi \hfill \end{array}\right]\text{.}$

For obtaining the stability condition, we have to find out eigenvalues of the amplification matrix T, but this study is very difficult since all the elements of T are different. Hence, the problem requires that the Eckert number E_c is assumed to be very small, that is, tends to zero. Under this consideration, we have E = 0, and the amplification matrix becomes

$$T=\left[\begin{array}{c}\hfill CD0\hfill \\ \hfill 0B0\hfill \\ \hfill 0KJ\hfill \end{array},\right]\text{.}$$

After simplification of the matrix T, we get the following eigenvalues:λ₁ = C, λ₂ = B, and λ₃ = J.For stability, each eigenvalue (λ₁, λ₂, and λ₃) must not exceed unity in modulus. Hence, the stability condition is|C| ≤ 1, |B| ≤ 1, and J ≤ 1, for all a, β.

Now, we assume that U is everywhere non-negative and V is everywhere non-positive. Thus, $B=1-2\left[a+b+2c\left(\frac{1+R}{P_{\text{r}}}+N_{\text{b}}\bar{C}+N_{\text{t}}\bar{T}\right)\right]\text{,}$ where $a=U\frac{\text{Δ}\tau }{\text{Δ}X},b=\left|V\right|\frac{\text{Δ}\tau }{\text{Δ}Y}\text{and}c=\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\text{.}$

The coefficients a, b, and c are all real and non-negative. We can demonstrate that the maximum modulus of B occurs when α ΔX = mπ and β ΔY = nπ, where m and n are integers; hence, B is real. The value of |B| is greater when both m and n are odd integers.

To satisfy the second condition |B| ≤ 1, the most negative allowable value is B = - 1. Therefore, the first stability condition is

$$2\left[a+b+2c\left(\frac{1+R}{P_{\text{r}}}+N_{\text{b}}\bar{C}+N_{\text{t}}\bar{T}\right)\right]\le 2\text{,}$$

(32)

that is,

$$\begin{array}{c}U\frac{\text{Δ}\tau }{\text{Δ}X}+\left|V\right|\frac{\text{Δ}\tau }{\text{Δ}Y}+\frac{2\left(1+R\right)}{P_{\text{r}}}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\hfill \\ +2N_{\text{b}}\bar{C}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}+2N_{\text{t}}\bar{T}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\le 1.\hfill \end{array}$$

(33)

Likewise, the third condition |J| ≤ 1 requires that

$$U\frac{\text{Δ}\tau }{\text{Δ}X}+\left|V\right|\frac{\text{Δ}\tau }{\text{Δ}Y}+\frac{2}{L_{\text{e}}}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\le 1.$$

(34)

Therefore, the stability conditions of the method are

$$\begin{array}{c}U\frac{\text{Δ}\tau }{\text{Δ}X}+\left|V\right|\frac{\text{Δ}\tau }{\text{Δ}Y}+\frac{2\left(1+R\right)}{P_{\text{r}}}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}+2N_{\text{b}}\bar{C}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\hfill \\ +2N_{\text{t}}\bar{T}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\le 1\text{and}U\frac{\text{Δ}\tau }{\text{Δ}X}+\left|V\right|\frac{\text{Δ}\tau }{\text{Δ}Y}\hfill \\ +\frac{2}{L_{\text{e}}}\frac{\text{Δ}\tau }{{\left(\text{Δ}Y\right)}^2}\le 1.\hfill \end{array}$$

Since, from the initial condition, $U=V=\bar{T}=\bar{C}=0$ at τ = 0, the consideration due to stability and convergence analysis is E_c < < 1 and R ≥ 0.5. Hence, convergence criteria of the method are P_r ≥ 0.37 and L_e ≥ 0.25.

Results and discussion

In order to investigate the problem under consideration, the results of numerical values of non-dimensional velocity, temperature, and species concentration within the boundary layer have been computed for different values of magnetic parameter M, radiation parameter R, Prandtl number P_r, Eckert number E_c, Lewis number L_e, Brownian motion parameter N_b, thermophoresis parameter N_t, and Grashof number G_r, respectively. To obtain the steady-state solutions of the computation, the calculations have been carried out up to non-dimensional time τ = 5 to 80. The velocity, temperature, and concentration profiles do not show any change after non-dimensional time τ = 50. Therefore, the solution for τ ≥ 50 is the steady-state solution. The graphical representation of the problem has been shown in Figures 3, 4, 5, 6, 7, 8, 9, and 10.

Figure 3

Grashof number ( G _r ) effect on velocity profiles.

Figure 4

Magnetic parameter ( M ) effect on velocity profiles.

Figure 5

Radiation parameter ( R ) effect on temperature profiles.

Figure 6

Eckert number ( E _c ) effect on temperature profiles.

Figure 7

Prandtl number ( P _r ) effect on temperature profiles.

Figure 8

Brownian parameter ( N _b ) effect on temperature profiles.

Figure 9

Thermophoresis parameter ( N _t ) effect on concentration profiles.

Figure 10

Lewis number ( L _e ) effect on concentration profiles.

In order to assess the accuracy of the numerical results, the present results (non-similar solution) are compared with the result obtained by Khan and Pop [19] (similar solution) and the values of magnetic parameter M, radiation parameter R, Eckert number E_c, and Grashof number G_r are considered zero (Table 1). From the comparison, excellent agreement is observed.

Table 1 Comparison of results for the reduced Nusselt number $N_{\text{u}}=-\frac{1}{\text{Δ}T}{\left(\frac{\partial T}{\partial y}\right)}_{y=0}$ when M  =  G _r  =  R  =  E _c  = 0 and P _r  =  L _e  = 10

In Figures 3, 4, 5, 6, 7, 8, 9, and 10, the dimensionless velocity, temperature, and concentration distributions are plotted against Y for the different non-dimensional time τ = 5 to 50 and corresponding values of Grashof number G_r, magnetic parameter M, radiation parameter R, Eckert number E_c, Prandtl number P_r, Brownian motion parameter N_b, thermophoresis parameter N_t, and Lewis number L_e, respectively.

In Figures 3 and 4, the velocity distribution is plotted respectively for different values of G_r and M. The non-dimensional time considered here is τ = 5, 20 and 50 and displays the entire step with a different pattern. Here, it is observed that when the values of G_r increase, then the velocity profiles increase and when the values of M increase, then the velocity profiles decrease.

In Figures 5, 6, 7, and 8, the temperature distribution is plotted respectively for different values of R, E_c, P_r, and N_b. The non-dimensional time considered here is τ = 5, 20 and 50 and displays the entire step with a different pattern. Here, it is observed that when the values of R increase, then the temperature profiles increases; when the values of E_c increase, then the temperature profiles also increase; when the values of P_r increase, then the temperature profiles decrease; and when the values of N_b increase, then the temperature profiles increase.

In Figures 9 and 10, the concentration distribution is plotted respectively for different values of N_t and L_e. The non-dimensional time considered here is τ = 5, 20 and 50 and displays the entire step with a different pattern. Here, it is observed that when the values of N_t increase, then the concentration profiles increase, but when the values of L_e increase, then the concentration profiles decrease.

From the present results and the result obtained by Khan and Pop [19], it was observed that the flow field shows the same trend with the variation of magnetic parameter M, radiation parameter R, Prandtl number P_r, Eckert number E_c, Lewis number L_e, Brownian motion parameter N_b, thermophoresis parameter N_t, and Grashof number G_r. However, the important part of this work is its comparison with the previous work, i.e., the present study is the unsteady case of Khan and Pop's [19] study when the values of magnetic parameter M, radiation parameter R, Eckert number E_c, and Grashof number G_r are considered zero.

Conclusions

An unsteady free convection boundary-layer flow of a nanofluid due to a stretching sheet is studied with the influence of magnetic field and thermal radiation. The explicit finite difference [24] technique with stability and convergence analysis has been employed as a solution technique to complete the formulation of the unsteady model. For the unsteady case (time-dependent), the non-dimensional time considered here is τ = 5, 20 and 50 and displays with the entire step and a different line pattern. The results are presented for the effect of various parameters. The velocity, temperature, and concentration effects on the sheet are studied and shown graphically.

From the present study, the concluding remarks have been taken as follows:

• 1. Larger values of the Grashof number showed a significant effect on momentum boundary layer.

• 2. The effect of the Brownian motion and thermophoresis stabilizes the boundary layer growth.

• 3. The boundary layers are highly influenced by the Prandtl number.

• 4. Using magnetic field, the flow characteristics could be controlled.

• 5. The thermal boundary layer thickness increases as a result of increasing radiation.

• 6. The presence of heavier species (large Lewis number) decreases the concentration in the boundary layer.

• 7. The Eckert number has a significant effect on the boundary layer growth.

Authors’ information

AI and LEA are assistant professors of Mathematics. MSK and IK are researchers.