Optimal and Numerical Solutions for an MHD Micropolar Nanofluid between Rotating Horizontal Parallel Plates.

The present analysis deals with flow and heat transfer aspects of a micropolar nanofluid between two horizontal parallel plates in a rotating system. The governing partial differential equations for momentum, energy, micro rotation and nano-particles concentration are presented. Similarity transformations are utilized to convert the system of partial differential equations into system of ordinary differential equations. The reduced equations are solved analytically with the help of optimal homotopy analysis method (OHAM). Analytical solutions for velocity, temperature, micro-rotation and concentration profiles are expressed graphically against various emerging physical parameters. Physical quantities of interest such as skin friction co-efficient, local heat and local mass fluxes are also computed both analytically and numerically through mid-point integration scheme. It is found that both the solutions are in excellent agreement. Local skin friction coefficient is found to be higher for the case of strong concentration i.e. n=0, as compared to the case of weak concentration n=0.50. Influence of strong and weak concentration on Nusselt and Sherwood number appear to be similar in a quantitative sense.

Introduction

The idea of micropolar fluid was introduced by Eringen [1-2]. This idea is a substantial generalization of the classical Navier Stokes model to discuss certain complex fluids. This particular class of fluids consists of rigid, randomly oriented spherical particles with microstructures such as liquid crystals, colloidal fluids, polymeric suspensions, hematological suspensions and animal blood etc. Extensive uses of micropolar fluid theory are given in the books of Lukaszewicz [3] and Eringen [4]. Bhargava et al [5] presented finite element solutions for mixed convective micropolar flow driven by a porous stretching sheet. Later on, Takhar et al [6] highlighted free convection MHD micropolar fluid flow between two porous vertical plates and perceived that velocity decreases with an increase in Hartman number. Stagnation point flow of a micropolar fluid towards a stretching sheet was investigated by Nazar et al [7]. Ziabakhsh et al [8] presented Homotopy analysis solutions of micropolar flow in a porous channel with heat and mass transfer. Similarly, Ishak et al [9] considered magnetohydrodynamic flow of micropolar fluid towards a stagnation point on a vertical surface. Joneidi et al [10] inspected behavior of micropolar flow in a porous channel with high mass transfer. Influence of chemical reaction and thermal radiation on MHD micropolar flow over a vertical moving porous plate in a porous medium with heat generation was discussed by Mohamed et al [11] who found that the translational velocity across the boundary layer and the magnitude of micro rotation at the wall decreased with an increase in magnetic field and Prandtl number. Some noteworthy studies on micropolar fluids with certain physical constraints can be found in [12-16].

In the last few decades nanofluids have proved to be extremely promising heat transfer agents in modern day industry and in numerous engineering applications of global interest like biomedical, optical and electronic fields. Choi [17] initiated the astonishing idea of nanofluids. Later numerous researchers and scientists investigated various real life flow problems under the influence of nanofluids. Buongiorno [18] presented a novel study on describing flow and heat transfer mechanisms of nanofluids. Similarly Nield et al [19] discussed natural convective boundary-layer flow of a nanofluid past a vertical plate and observed that the reduced Nusselt number is a decreasing function of thermophoresis and Brownian motion parameter. Nadeem et al [20] presented the optimized analytical solution for oblique flow of a Casson nanofluid with convective boundary conditions. They found that nanoparticle concentration is an increasing function of stretching parameter and Brownian motion while it is a decreasing function of thermophoresis, Biot number and non-Newtonian (Casson) parameter. Ganji et al [21] conducted a valuable study on heat transfer of Cu-Water nanofluid flow between parallel plates and concluded that heat flux at the surface has a direct relationship with nanoparticle volume fraction. Ganji et al [22] discussed simulation of MHD Cuo-Water nanofluid flow and convective heat transfer under the influence of Lorentz forces. Some notable studies on the topic can be found in [23-41].

The main goal of the present study is to discuss hydromagnetic flow of a micropolar nanofluid between parallel plates. Optimal homotopy analysis method proposed by Liao [24-25] is utilized to obtain graphical results against embedding physical parameters. Numerical values for skin friction, Nusselt and Sherwood numbers are obtained through OHAM and also numerically through mid-point integration scheme [26-27]. It would be obligatory to mention that mid-point integration scheme is used to verify our obtained analytical results. It is found that both the results are in a very good agreement with each other.

Mathematical Formulation

Consider steady incompressible 3D flow of an electrically conducting micropolar nanofluid between two horizontal parallel plates. Both the fluid and the plates rotate together around the y-axis with a constant angular velocity Ω. The x and y − axes are taken along and perpendicular to the plate respectively while the z − axis is taken normal to the xy plane. The plates are placed at y = 0 and y = h. The lower plate is being stretched by two equal and opposite forces so that the position of the point (0,0,0) remains unchanged. A uniform magnetic flux with density B 0 acts along the y-axis. The upper plate is subject to constant wall suction velocity v 0(< 0) or constant wall injection velocity v 0(> 0) respectively as shown in Fig 1.

Fig 1

Physical description of the problem.

The governing equations of the flow problem can be stated as [5–11, 18]

The appropriate boundary conditions are where u, v and w are the velocity components along the x, y, z — directions respectively while u = ax shows that lower plate is stretching.ρ is the density, ϑ is the kinematic viscosity, is the vortex viscosity, N is the micro rotation velocity, Ω is the rotation velocity, B 0 is the magnetic field, σ is electrical conductivity, α is thermal diffusivity, k is the thermal conductivity, C is the specific heat of the fluid, j is the micro-inertia density, v 0 is the suction/injection velocity, T is the temperature of the fluid, T 0 is the temperature at lower plate, T is the temperature at upper plate. D Brownian diffusion coefficient, D is the thermophoresis diffusion coefficient, C is the concentration and n is the boundary parameter. The case n = 0 represents strong concentration, n = 0.5 indicates weak concentration and vanishing of the antisymmetric part of the stress tensor and n = 1 represents turbulent flow. Spin gradient viscosity γ is defined as

Using transformations

Eq (1) is identically satisfied while Eqs. (2)–(7), in dimensionless form are given by with boundary conditions where primes denote differentiation with respect to η. N 1 is the coupling parameter, N 2 is the spin-gradient viscosity parameter, R is the Reynold number, M is the magnetic parameter, Kr is the rotation parameter, Pr is the Prandtl number, Sc is the Schmidt number, N is the Brownian motion parameter and N is the thermophoresis parameter which are defined by

The skin friction coefficient C , Nusselt number Nu and Sherwood number Sh are defined as where

The dimensionless forms of skin friction coefficient Nusselt number Nu and Sherwood number Sh are

Method of Solution

1. Optimal HAM Solution

The governing system of coupled ordinary differential Eqs (10)–(14) is nonlinear and extremely complicated in nature. These equations are solved analytically by optimal homotopy analysis method (OHAM). Following Liao [24-25] we know that f(η), g(η), θ(η), G(η) and ϕ(η) can be expressed by a set of a certain type of exponential base functions in the forms in which and are series coefficients. The initial guesses f 0, g 0, θ 0, G 0 and ϕ 0 for f(η), g(η), θ(η), G(η) and ϕ(η) are selected as follows

We select the auxiliary linear operators as

The above linear operators have the following properties Where C 0 − C 11 are arbitrary constants.

The zeroth order homotopic deformation equations can be written as where p ϵ [0, 1] indicates the embedding parameter and and the nonzero auxiliary parameters. Moreover the nonlinear operators N , N , N , N and N are prescribed as

When p varies from 0 to 1, we have

By means of Taylor's series

The auxiliary converging parameters are chosen in such a way that the series (50) converges when p = 1. Thus we have

The resulting mth – order deformation equations are with the following definitions in which

The general solutions of Eqs (56)–(60) can be written as in which and are the particular solutions of the Eqs (68)–(72) Note that Eqs (56)–(60) can be solved by Mathematica one after the other in the order m = 1,2,3, …

Optimal convergence-control parameters

We know that the homotopic series solutions involve the non-zero auxiliary parameters and which define the convergence-region of the Homotopy series solutions. To determine the optimal values of and we utilize the concept of minimization by defining the average squared residual errors as introduced by [24].

Following Liao [24]

where is total averaged squared residual error, δη = 0.5, k = 20. Total and individual averaged squared residual errors are computed using N 1 = N 2 = N 3 = 0.10, R = 0.20, n = 0.50, M = 0.10 = Kr = Nt = Nb = λ, Pr = 1 = Sc. By means of computational software Mathematica we obtained the total and individual average squared residual errors at various order of iterations using highly efficient Mathematica package . 0 which can be found at http://numericaltank.sjtu.edu.cn/BVPh2_0.htm. The basic idea is to minimize the total average squared residuals and determining the corresponding local optimal convergence control parameters. For this purpose, Tables 1 and 2 are prepared for the case of various optimal convergence control parameters. Table 1 gives the minimum value of total averaged squared residual error at several order of iterations while Table 2 is presented for the individual average squared residual error at different orders of approximations using the optimal values from Table 1 at m = 6. It is quite obvious from these two tables that the averaged squared residual errors and total averaged squared residual errors continue decreasing as we increase the order of approximation. So, Optimal Homotopy Analysis Method is an extremely effective tool in obtaining convergent series solutions for highly nonlinear systems of differential equations.

Table 1

Optimal convergence control parameters and total averaged squared residual errors using BVPh2. 0 when N1 = N2 = N3 = 0.10, R = 0.20, n = 0.50, M = 0.10 = Kr = Nt = Nb = λ, Pr = 1 = Sc.

m	c0f	c0g	c0θ	c0G	c0ϕ	Emt	CPU TIME[S]
2.0	−0.79	−0.74	−1.06	−7.47	−0.83	5.86 × 10−3	3.88
4.0	−0.87	−0.76	−1.10	−8.16	−0.89	4.78 × 10−6	116.0
6.0	−0.89	−0.88	−0.91	−8.28	−1.08	1.67 × 10−9	2453

Table 2

Individual Averaged squared residual errors using optimal values at m = 6 from Table 1.

m	Emf	Emg	Emθ	EmG	Emϕ	CPUTIME[S]
4.0	1.04 × 10−5	1.17 × 10−8	3.50 × 10−9	6.18 × 10−8	1.30 × 10−7	2.12
6.0	1.59 × 10−9	1.05 × 10−12	1.13 × 10−12	3.80 × 10−11	3.96 × 10−11	3.98
10	7.40 × 10−16	2.38 × 10−20	1.46 × 10−19	2.39 × 10−17	4.57 × 10−18	8.69
14	4.85 × 10−22	3.44 × 10−28	1.91 × 10−26	1.58 × 10−23	4.31 × 10−25	15.28

2. Numerical Solution

Despite OHAM, the governing system of Eqs (10)–(14) is also solved numerically using midpoint integration as a basic scheme and Richardson extrapolation as an enhancement scheme with highly efficient computational software Maple as used by several other authors [26-30] This scheme works by transforming the governing system of nonlinear higher order differential equations into a system of first order ordinary differential equations which are then solved using iterative schemes through midpoint integration as it belongs to the class of higher order Runge-Kutta methods. A mesh size of Δh = 0.001 was set for a convergence criterion of 10−6 in our computations. Our computed numerical results are in very good agreement with analytical results obtained by optimal homotopy analysis method. We are thus confident that our applied numerical algorithm is up to the mark.

Results and Discussion

The aim here is to discuss the behavior of velocity, temperature, micro rotation and concentration profiles against emerging physical parameters in our flow problem. Figs 2–25 are plotted for this purpose. Figs 2–7 are plotted to determine the influence of coupling parameter N 1, spin gradient viscosity parameter N 2, rotation parameter Kr, Hartman number M, porosity parameter λ and Reynolds number R on velocity profile f′(η). From Figs 2 and 3 we observe that the influence of coupling parameter N 1 and spin gradient viscosity parameter N 2 on the velocity profile f′(η) is similar, i.e. it initially decreases near the lower stretching plate and from the center of the plates towards the upper plate the behavior is reversed.

Fig 2

Effect of N 1 on f′(η).

Fig 25

Effect of Sc on ϕ(η).

Fig 7

Effect of R on f′(η).

Fig 3

Effect of N 2 on f′(η).

On the other hand, from Figs 4 and 5 velocity profile f′(η) increases towards the main stream with an increase in rotation parameter Kr Hartman number M and then it tends to decrease towards the upper plate.

Fig 4

Effect of Kr on f′(η).

Fig 5

Effect of M on f′(η).

The effects of porosity parameter λ and Reynolds number R are revealed through Figs 6 and 7. From Fig 6 we observe that increasing the porosity parameter λ results in an increase in the velocity profile f′(η) from lower towards the upper plate. From Fig 7 velocity profile f′(η) decreases near the lower plate but increases towards the upper plate.

Fig 6

Effect of λ on f′(η).

The impact of embedding physical parameters such as coupling parameter N 1, rotation parameter Kr, Hartman number M, porosity parameter λ and Reynolds number R on transverse velocity g(η) is displayed through Figs 8–12. Figs 8 and 9 depict that coupling parameter N 1 and rotation parameter Kr have opposite influence on transverse velocity g(η). Similarly from Fig 10 increasing the magnetic field through Hartman number M increases the transverse velocity g(η) near the lower plate but contrary behavior is observed near the upper plate. It is quite evident from Figs 11 and 12 that increasing the porosity parameter λ and Reynolds number R have a positive influence on the transverse velocity g(η).

Fig 8

Effect of N 1 on g(η).

Fig 12

Effect of R on g(η).

Fig 9

Effect of Kr on g(η).

Fig 10

Effect of M on g(η).

Fig 11

Effect of λ on g(η).

Variation in the micro rotation profile G(η) against various parameters such as coupling parameter N 1, spin gradient viscosity parameter N 2, Hartman number M and Reynolds number R are presented through Figs 13–16. From Figs 13 and 14 the effects of coupling parameter N 1 and spin gradient viscosity parameter N 2 on the micro rotation profile G(η) are quite opposite to each other. From Fig 13 increasing the coupling parameter N 1 leads to an increase in the micro rotation profile G(η) while it decreases with an increase in spin gradient viscosity parameter N 2 (See Fig 14). Influence of Hartman number M and Reynolds number R on the micro rotation profile G(η) is found to be similar from Figs 15 and 16. Both these parameters enhance the micro ration profile G(η) between the parallel plates.

Fig 13

Effect of N 1 on G(η).

Fig 16

Effect of R on G(η).

Fig 14

Effect of N 2 on G(η).

Fig 15

Effect of M on G(η).

The impact of physical parameters such as Brownian motion parameter N , thermophoresis parameter N , Prandtl number Pr and Reynolds number R on the temperature profile θ(η) is presented through Figs 17–20. From Figs 17 and 18 temperature profile θ(η) increases significantly with an increase in Brownian motion N as well as thermophoresis parameter N . Figs 19 and 20 illustrate effects of Prandtl number Pr and Reynolds number R on the temperature profile θ(η). Temperature profile θ(η) increases with an increase in Prandtl number Pr(Fig 19), while it decreases with increasing Reynolds number R (Fig 20).

Fig 17

Effect of N on θ(η).

Fig 20

Effect of R on θ(η).

Fig 18

Effect of N on θ(η).

Fig 19

Effect of Pr on θ(η).

The effect of physical parameters such as Brownian motion parameter N , thermophoresis parameter N , Prandtl number Pr, Reynolds number R and Schmidt number Sc on the concentration profile φ(η) is presented through Figs 21–25. Figs 21 and 22 indicate that Brownian motion parameter N increases while thermophoresis parameter N decreases the concentration profile φ(η) between the plates. From Figs 23 and 24 concentration profile φ(η) decreases with an increase in Prandtl number Pr while it increases with an increase in Reynolds number R. Finally Fig 25 shows that for large values of Schmidt number Sc, concentration profile φ(η) decreases between the parallel plates.

Fig 21

Effect of N on ϕ(η).

Fig 22

Effect of N on ϕ(η).

Fig 23

Effect of Pr on ϕ(η).

Fig 24

Effect of R on ϕ(η).

Physical quantities of interest such as skin friction co-efficient, Nusselt and Sherwood number are computed through Tables 3–5. Table 3 examines the influence of embedding parameters such as N 1, N 3, M, Kr and R on skin friction coefficient. All these parameters have a positive influence on skin friction at the lower plate but it is interesting to mention here that this increasing behavior is more prominent for the case of strong concentration i.e. n = 0.0 when compared to the case of weak concentration n = 0.5. The impact of coupling parameter N 1, Brownian motion parameter N , thermophoresis parameter N , Reynolds number R and Schmidt number Sc on the Nusselt number is depicted through Table 4. For higher values of coupling parameter N 1, Nusselt number –θ′(0) tends to increase while it decreases with an increase in, Brownian motion parameter N , thermophoresis parameter N . Similarly, it is also observed that with an increase in Reynolds number R, Nusselt number –θ′(0) also increases while it decreases with an increase in Schmidt number Sc for strong n = 0.0 as well as weak concentration n = 0.5. Finally Table 5 depicts that Sherwood number –ϕ′(0) responds in an opposite manner against Brownian motion parameter N , thermophoresis parameter N . Moreover local mass flux –ϕ′(0) decreases with Reynolds number R and it increases with an increase in Schmidt number Sc for strong n = 0.0 as well as weak concentration n = 0.5.

Table 3

Numerical values of skin friction at the wall when N 2 = 1.0.

					Cf = (1 + N 1)f "(0) + N 1 G(0)
					n = 0.0 (Strong Concentration)		n = 0.50 (Weak Concentration)
N 1	N 3	M	Kr	R	Numerical	Optimal Ham	Numerical	Optimal Ham
0.0	0.1	0.1	0.1	0.5	−3.46067	−3.46067	−3.46067	−3.46067
0.1					−3.79947	−3.79947	−3.62896	−3.62896
0.5					−5.13213	−5.13213	−4.31894	−4.31894
0.1	0.0	0.1	0.1	0.5	−3.79947	−3.79947	−3.62881	−3.62881
	0.3				−3.79947	−3.79947	−3.62924	−3.62924
	0.5				−3.79947	−3.79947	−3.62953	−3.62953
0.1	0.1	0.0	0.1	0.5	−3.78723	−3.78723	−3.61727	−3.61727
		0.5			−3.84813	−3.84813	−3.67540	−3.67540
		1.0			−3.90829	−3.90829	−3.73284	−3.73284
0.1	0.1	0.1	0.0	0.5	−3.79942	−3.79942	−3.62890	−3.62890
			0.5		−3.80078	−3.80078	−3.63020	−3.63020
			1.0		−3.80486	−3.80486	−3.63410	−3.63410
0.1	0.1	0.1	0.1	0.5	−3.79947	−3.79947	−3.62896	−3.62896
				1.0	−3.84769	−3.84769	−3.67514	−3.67514
				1.5	−3.89581	−3.89581	−3.72123	−3.72123

Table 5

Numerical values of Mass flux at the wall when N1 = N2 = 1.0, N3 = 0.1, λ = 0.1, M = 0.1, Kr = 0.1, Pr = 1.2.

				Sh = −φ ′(0)
				n = 0.0 (Strong Concentration)		n = 0.50 (Weak Concentration)
N b	N t	R	Sc	Numerical	Optimal Ham	Numerical	Optimal Ham
0.5	0.1	0.5	1.0	1.09115	1.09115	1.09058	1.09058
0.8				1.08604	1.08604	1.08571	1.08571
1.2				1.08064	1.08064	1.08030	1.08030
0.1	0.2	0.5	1.0	1.30437	1.30437	1.30476	1.30476
	0.5			2.47999	2.47999	2.48122	2.48122
	0.8			4.41360	4.41360	4.41532	4.41532
0.1	0.1	0.5	1.0	1.11270	1.11270	1.11274	1.11274
		1.0		1.10926	1.10926	1.10935	1.10935
		1.5		1.10493	1.10493	1.10509	1.10509
0.1	0.1	0.5	0.5	1.09558	1.09558	1.09580	1.09580
			1.0	1.11270	1.11270	1.11274	1.11274
			1.5	1.12995	1.12995	1.12980	1.12980

Table 4

Numerical values of heat flux at the wall when N2 = 1.0, N3 = 0.1, λ = 0.1, M = 0.1, Kr = 0.1, Pr = 1.

					Nu = −θ ′(0)
					n = 0.0 (Strong Concentration)		n = 0.50 (Weak Concentration)
N 1	N b	N t	R	Sc	Numerical	Optimal Ham	Numerical	Optimal Ham
0.0	0.1	0.1	0.5	1.0	0.92042	0.92042	0.92042	0.92042
0.1					0.92044	0.92044	0.92043	0.92043
0.3					0.92049	0.92049	0.92043	0.92043
0.1	0.5	0.1	0.5	1.0	0.71075	0.71075	0.71074	0.71074
	0.8				0.57831	0.57831	0.57831	0.57831
	1.2				0.43240	0.43240	0.43240	0.43240
0.1	0.1	0.5	0.5	1.0	0.71285	0.71285	0.71284	0.71284
		0.8			0.58125	0.58125	0.58124	0.58124
		1.2			0.43578	0.43578	0.43575	0.43575
0.1	0.1	0.1	0.5	1.0	0.92044	0.92044	0.92043	0.92043
			1.0		0.95655	0.95655	0.95653	0.95653
			1.5		0.99309	0.99309	0.99306	0.99306
0.1	0.1	0.1	0.5	0.5	0.92078	0.92078	0.92077	0.92077
				1.0	0.92044	0.92044	0.92043	0.92043
				1.5	0.92010	0.92010	0.92009	0.92009

Concluding Remarks

Hydromagnetic flow of micropolar nanofluid between two horizontal parallel plates in a rotating system has been investigated numerically as well as analytically using optimal HAM. The core outcomes of the study can be enlisted as

The effects of coupling parameter N 1 and Hartman number M on velocity profile f′(η) are found to be opposite compared to viscosity parameter N 2 and rotation parameter Kr.

Influence of magnetic field M and coupling parameter N 1 on transverse velocity g(η) are opposing each other. Moreover it is observed that increasing the porosity parameter λ leads to an increase in the transverse velocity g(η).

It is observed that temperature profile θ(η) enhances for large values of Prandtl number Pr, thermophoresis parameter Nt and Brownian motion parameter Nb for the case of strong concentration (n = 0).

Micro rotation profile G(η) raises with coupling parameter N 1, magnetic field parameter M and Reynolds number R for strong concentration i.e. when n = 0.

Increasing the thermophoresis parameter Nt and Brownian motion parameter Nb results in a contrary behavior on concentration profile ϕ(η).

Skin friction at the lower plate tends to increase with magnetic field parameter M, rotation parameter Kr and Reynolds number R. This behavior is found to be more prominent for the case of strong concentration i.e. (n = 0) compared to the weak concentration (n = 0.5).

Heat flux –θ′(0) at the lower plate drops with Increasing thermophoresis parameter Nt and Brownian motion parameter Nb while on the other hand mass flux –ϕ′(0) responds in an opposite manner against these two parameters.