Numerical Solution of Nonlinear Diff. Equations for Heat Transfer in Micropolar Fluids over a Stretching Domain

A numerical study based on finite difference approximation is attempted to analyze the bulk flow, micro spin flow and heat transfer phenomenon for micropolar fluids dynamics through Darcy porous medium. The fluid flow mechanism is considered over a moving permeable sheet. The heat transfer is associated with two different sets of boundary conditions, the isothermal wall and isoflux boundary. On the basis of porosity of medium, similarity functions are utilized to avail a set of ordinary differential equations. The non-linear coupled ODE’s have been solved with a very stable and reliable numerical scheme that involves Simpson’s Rule and Successive over Relaxation method. The accuracy of the results is improved by making iterations on three different grid sizes and higher order accuracy in the results is achieved by Richardson extrapolation. This study provides realistic and differentiated results with due considerations of micropolar fluid theory. The micropolar material parameters demonstrated reduction in the bulk fluid speed, thermal distribution and skin friction coefficient but increase in local heat transfer rate and couple stress. The spin behavior of microstructures is also exhibited through microrotation vector N(η).


Introduction
The scientific and technological advances have brought a great awakening of interest in constructing different types of fluids and investigating their flow behavior in various useful geometries. Fluids with microstructures behave differently from the classical fluids. The flow and heat transfer behavior of these fluids cannot be described adequately with the classical theory of Newtonian fluid flows. Several theories have been presented to describe the very nature of these fluids. However, theory of micropolar fluids presented by Eringen [1] provides ample details required for justification of dynamics of such fluids. Micropolar fluids consist of rigid, randomly oriented, spherical particles with their own spins and microrotation, suspended in a viscous medium. Here the microelements are allowed to undergo rigid rotations only without stretch. The micropolar fluid model, apart from usual velocity vector involves a microrotation vector and a gyration parameter to simulate the kinematics of microrotation. These fluids possess monosymmetric stress tensor. Later, Eringen [2] extended his theory for thermo-micropolar fluids and derived the constitutive laws. Ariman et al. [3,4] presented an excellent review of micropolar fluids and their applications. Ahmadi [5] investigated the boundary layer flow of micropolar fluids past a semi-infinite plate. The basic theory of micropolar fluids can be viewed in the book written by Eringen [6] as well as by Be'g et al. [7]. Rehman et al. [8] considered heat transfer in two-dimensional steady hydromagnetic natural convection flow of a micropolar fluid past a non-linear stretching sheet with temperature dependent viscosity.
A porous medium is usually composed of a solid matrix and voids, it has relevance to heat and mass transfer applications in thermal and geophysical processes. The viscous flow through porous medium is described by the well-known Darcy law which was supported experimentally, numerically and theoretically [9,10] but as stated above, it is valid only for flows with a sufficiently small velocity or R Reynolds number. The Darcy law has been formulated ∇p = − µ K u by Bear [11]. The constant in Darcy's equation was proved later by Muskat and Wyckoff [12] that is related to the permeability of the porous material. Liu [13] presented closed form solution for the flow and heat transfer of a viscous fluid saturated in a porous medium past a permeable and non-isothermal stretching sheet with internal heat generation or absorption and radiation. Mohamed and Kasseb [14] obtained numerical solution steady flow and heat transfer in a porous medium saturated with a Sisko nano fluid (non-Newtonian power-law) over a nonlinearly stretching sheet in the presence of heat generation/absorption. Ferdows et al. [15] investigated effects of thermal radiation on a steady boundary layer flow with temperature-dependent thermal conductivity due to a stretching sheet through porous medium the presence of transverse magnetic field near a stagnation point. Berre et al. [16] described detailed reviews on modeling approaches for flow in fractured porous media, from physical, conceptual and mathematical models with two discretization approaches. Sun [17] presented an extension of the explicit moving particle semi-implicit (MPS) method with model formulation, based on the local volume averaging equations to compute the macroscopic behaviors of incompressible flows in porous media.
Heat transfer is important in many industrial processes. Upendar and Srinivasacharya [18] analyzed a mathematical model for the steady, mixed convection heat and mass transfer along a semi-infinite vertical plate embedded in a micropolar fluid in the presence of a first-order chemical reaction and radiation. MHD flow of micropolar fluid past a stretching sheet with heat transfer and with suction/blowing through a porous medium has been studied by Aldabe et al. [19]. Sharma et al. [20] studied the fully developed electrically conducting micropolar fluid flow and heat transfer along a semi-infinite vertical porous moving plate including the effect of viscous heating and in the presence of a magnetic field applied transversely to the direction of the flow. Mohammedein and Gorla [21] analyzed the flow of micropolar fluids bounded by a stretching sheet with a prescribed wall heat flux, viscous dissipation and internal heat generation. Abo-Eldahab and El-Aziz [22] considered heat transfer effect in a micropolar fluid flow induced by a stretching surface immersed in a porous medium with uniform free stream. Ahmad et al. [23] obtained closed form solution for a viscous, incompressible, MHD flow over a porous stretching sheet. Mahapatra and Gupta [21] investigated the flow and heat transfer characteristics over a stretching sheet with a uniform magnetic field and prescribed surface heat flux. Dayyan et al. [24] studied the Newtonian fluid flow with heat transfer through porous medium and presented analytical solution by employing the Homotopy Analysis Method (HAM).
The mathematical formulation of fluid flow problems is based on the fundamental laws of conservation and the nonlinear characteristic is inherited in the governing equations. Thus, in general, the exact solution for the flow problems is difficult to obtain, as opined by Ren [25]. However, different numerical and analytical approximate methods have been adopted to analyze the nonlinearity in the flow process. Some most often used approximate analytical methods for analysis of flow problems include VIM (variational iteration method), HAM (homotopy analysis method), DTM (differential transformation method), ADM (Adomain decomposition method), VPM (variation of parameter method), OHAM (optimal homotopy asymptotic method), etc. However, the computational cost and time of these methods is increased for the determination of the unknowns (included to meet the second boundary conditions). Runge Kutta (R-K) methods with or without shooting techniques have been widely utilized for the solution of flow problems [26][27][28]. However, the Runge-Kutta method is widely used to solve the ordinary differential equation [29]. Similarly, other basic numerical methods have their own limitations. In the meantime, the advent of powerful computers facilitated the implementation of various numerical approaches to handle the complexities in this research area. In this scenario, numeric computation is capable of yielding reliable and cost-effective numerical solutions. Here, the governing differential equations are discretized. The basic discretization methods are FDM (finite difference method), FVM (finite volume method), and FEM (finite element method). FDM is the oldest of the three techniques, it is easy to implement and works well for the simple grids. In FDM, the continuous derivative is discretized to represent the difference of the flow variable between neighboring cells [30]. Among others, Ahmad et al. [31] used a finite-difference scheme known as the Keller-box method to study the flow and heat transfer of a micropolar fluid past a nonlinearly stretching plate. Finite difference solution for micropolar flow problems are analyzed by Hussain et al. [32] and Ashraf et al. [33], Shafique and Rashid [34].
Motivated from the above cited studies, we revisited [24] for the non-linear analysis of thermal transportation in micropolar fluids' motion over a permeable stretching sheet related to porous media with wall temperature and heat flux boundary conditions. A finite difference based scheme which comprises successive over relaxation iterative procedure, Simpson's (1/3) rule and Richardson extrapolation is harnessed to yield the solution. The optimum value of the relaxation parameter ω opt is estimated to accelerate the convergence of the SOR method. In computational terms, this scheme is cost effective and well established. It can be utilized for reliable solution of simultaneous non-linear equations. We obtained results for the velocity, skin friction coefficient temperature, local Nusselt number, microrotation and couple stress under the variation of controlling parameters. Section 2 presents mathematical analysis, Section 3 provides the detail of numerical method. Section 4 contains results with discussion. Accuracy for this solution is assured through their computations for three grids sizes (h, h/2, h/4) and their closeness to other published studies in the limiting case. In addition, the extrapolation routine raises the order of accuracy to higher level.

Mathematical Analysis
In the mathematical theory of micropolar fluids, there are, in general, six degrees of freedom, three for translation and three for microrotation. The essence of the theory of micropolarfluid flow lies in the extension of the constitutive equations for Newtonian fluids so that more complex with microstructure such as animal blood, muddy water, colloidal fluids, lubricants and chemical suspensions. The micropolar fluids consist of randomly oriented particles suspended in a viscous medium. In practice, the theory of micropolar fluids requires that one must add a transport equation representing the principle of conservation of local angular momentum to the usual transport equations for the conservation of mass and momentum, and additional local constitutive parameters are also introduced. In this way, it enables to recover the inadequacy of Navier-Stokes theory to describe the correct behavior of fluids with microstructures. The governing equations of the motion are [1]: where V is velocity, Ω is microrotation vectors, ρ is the density, f is body force, l is body couple per unit mass, respectively, p is pressure, j is the micro-inertia, µ, κ, λ 1 , γ are, respectively, dynamic viscosity, vortex viscosity, Stokes viscosity and the spin gradient viscosity. α , β , λ symbolize material constants.
Consider micropolar fluid flow through a homogeneous porous medium of permeability of K, over a porous stretching sheet. The flow is time independent, incompressible and two-dimensional. Fluid flows due to a permeable sheet of length L which is stretching linearly along the x − axis. The y − axis is perpendicular to the sheet. The linear velocity distribution of flow along the sheet is u w = u 0 x/L. The body force and body couple is neglected. The schematic diagram is shown in Figure 1.
where V is velocity, Ω is microrotation vectors, ρ is the density, f is body force, l is body couple per unit mass, respectively, p is pressure, j is the micro-inertia, μ , κ , 1 λ , γ are, respectively, dynamic viscosity, vortex viscosity, Stokes viscosity and the spin gradient viscosity.
Consider micropolar fluid flow through a homogeneous porous medium of permeability of K , over a porous stretching sheet. The flow is time independent, incompressible and two-dimensional.
Fluid flows due to a permeable sheet of length L which is stretching linearly along the axis x −  .
where v u, are velocity components along horizontal and vertical directions and Ω ) , is micro rotation vector perpendicular to − xy plane, having resemblance to ∇ × V . The fluid temperature is . T The energy equation: The hydrodynamic boundary conditions are: where 0 u is constant, w v is mass transfer speed at the wall and m is constant in the range The case m = 0 is called strong concentration by Guram and Smith [35], and indicates N = 0 near the surface and represents concentrated particle flows in which the microelements close to the surface are unable to rotate (Jena and Mathur [36]). The case m = 1/2 indicates the vanishing of antisymmetrical part of the stress tensor and denotes weak concentration (Ahmadi [5]). The case m = 1, as suggested by Peddieson [37], is used for the modeling of turbulent boundary-layer flows.
The thermal boundary conditions are: The Equations (1)-(3) are respectively simplified in view of the above assumptions to yield as below: where u, v are velocity components along horizontal and vertical directions and Ω= Ω(0, 0, ω 3 ) is micro rotation vector perpendicular to xy− plane, having resemblance to ∇ × V. The fluid temperature is T. The energy equation: The hydrodynamic boundary conditions are: where u 0 is constant, v w is mass transfer speed at the wall and m is constant in the range 0 ≤ m ≤ 1. The case m = 0 is called strong concentration by Guram and Smith [35], and indicates N = 0 near the surface and represents concentrated particle flows in which the microelements close to the surface are unable to rotate (Jena and Mathur [36]). The case m = 1/2 indicates the vanishing of antisymmetrical part of the stress tensor and denotes weak concentration (Ahmadi [5]). The case m = 1, as suggested by Peddieson [37], is used for the modeling of turbulent boundary-layer flows. The thermal boundary conditions are: where s is power law index and Tau (τ) is the effective thermal conductivity of the medium and is a function of thermal conductivities of the fluid and solid phases and the porous medium microstructure.
x * = x/L. Let us take the similarity transformations as: The Equation (4) is readily satisfied. The Equations (5)- (7) are respectively transformed to ordinary differential forms: (1 where prime denotes the differentiation with respect to η. R = ρu 0 K Lµ is the Reynolds number, ρjκµ γ are non-dimensional material constants. The boundary conditions (8) and (9) then become: Furthermore, λ > 0 shows suction and λ < 0 is for injection. The physical quantities of significance or local wall shear stress τ w , the wall couple stress m w , and the heat transfer from sheet surface q w as reported in [33].

Existence of Results
In order to establish the existence of a solution, the Equation (12) is formulated to the corresponding initial value problem (IVP), and the related boundary conditions in (15) are rewritten as: The parameter δ corresponds to coefficient of skin friction. This IVP formulation enables us to utilize topological shooting argument [38]. The solution of Equation (12) along with (16) is denoted by Then it is to show that a choice of δ can be made such that f (η, δ) exists, for all η > 0 and f (∞) = 0 is satisfied. To achieve this, two sets are defined: Two lemmas are presented below: The set X is nonempty and open. (12) when η = 0, to get
Next for X to be open, let δ ∈ X. We will show that all δ sufficiently close to δ are in X.
At η x (δ), we have 0 < f < 1 and f = 0. The Equation (12) at η x (δ) implies that Thus by continuity of the solution of the IVP in the initial condition, for δ sufficiently close to δ, f (η, δ) = 0 will also have a root near η x (δ) near η x (δ) with f (η, δ) > 0. Thus δ ∈ X. This leaves only the possibility that f = 0 and f = 0, simultaneously; however, setting this information to Equation It contradicts to Equation (16), where f = 1.

Lemma 2. The set Y is nonempty and open.
By integrating Equation (12) on [0, η], one is given Now it is to show that there are values of δ < 0 when magnitude of δ is large such that f = 0 in the interval (0, 1], say strictly before f = 0. Now suppose that this assertion is false and consider the following cases.
We now show that we must have C = 0; for this, we suppose that 0 <C< 1. Since f (η) < 0 for η > 0, f is bound below by C > 0 and so f tends to positive infinity. The term f f becomes negative. From Equation (12), we have There exists a point η 2 > 0 such that η > η 2 , we have . When η → ∞ implies that f → ∞ , contradicting the fact that f (η) < 0, so f (∞, δ * ) = 0, it establishes the theorem.

Numerical Method
The finally resulted set of non-linear coupled ordinary differential equations is solved numerically subject to physical boundary conditions. Firstly, the third order derivative is reduced to second order by setting: The Equation (12) is transformed as below: Suppose ϕ represents each of the dependent functions: f , q, N, θ and then using central differences defined below: The second order differential Equations (13), (14) and (18) are now respectively approximated by central difference at typical grid point η = η n of the interval [0, ∞). The discretization of the domain [0, ∞) is uniform with a step size h. We obtain The boundary conditions (15) become: The algebraic Equations (18)-(20) are solved iteratively by employing the SOR method (see Smith [39] p. 262) and the first order ODE (16) is integrated by Simpson's (1/3) rule Gerald [40] combined with a corrector formula (see Miline [41] p. 48) subject to the associated boundary conditions (21). The results for f (η), N(η), θ(η), − f (0), −N (0) and −θ (0) are in order of accuracy o(h 2 ) because of the inherent second order approximation of central finite differences. The accuracy of the scheme is checked with repetition of the procedure at the grid size h, h/2 and h/4. The results obtained at three grid sizes are utilized to involve Richardson Extrapolation to the limit (see Burden [42] page 168) to achieve accuracy in the results in the order o(h 6 ). The solution procedure, which is mainly based on the algorithm described in Syed et al. [43] is used to accelerate the iterative procedure and to improve the convergence. The relaxation parameter ω is optimized as proposed by Nakamura [44]. An initial guess is taken for ω as arbitrarily after each relaxation sweep except the first one, then the optimal values of ω are estimated with the formula given below: and where i and t, respectively, stand for the grid point and relaxation sweep. The process of optimization is stopped when the following criterion is met: subject to the conditions µ i ω < 1 and i > 1. The iterative process is continued with ω = ω opt and 10 −5 ≤ E opt ≤ 10 −3 . The iterative procedure is stopped if the following criteria is satisfied for four consecutive iterations here TOL iter is the prescribed error tolerance. The numerical scheme is coded in robust scientific programming language Fortron-90.

Results and Discussion
The theoretical and numerical analysis for this work was focused to reveal the physical nature of the micropolar fluids flow through porous medium caused by a permeable stretching sheet with both isothermal wall temperature condition and isoflux boundary condition. An extensive computational effort has been carried out to capture the reliable results for non-dimensional horizontal flow speed f (η), microrotation N(η) and temperature function θ(η), skin friction coefficient − f (0), heat transfer rate at surface −θ (0) and couple stress −N (0). The physical behavior of these quantities of interest has been scrutinized through the suitable ranges of the pertinent parameters-namely Reynolds number R, micropolar material parameter C 1 , suction/injection parameter λ, Prandtl number Pr, heat index parameter s. The non-dimensional micropolar material parameters C 1 , C 2 and C 3 are chosen arbitrarily. The values of C 1 , are considered as 0.5, 1.0, 1.5 and 4.0, while C 2 = 0.1 and C 3 = 0.5 are chosen arbitrarily. However, if C 1 = 0 and N = 0, this problem corresponds to Newtonian fluids flow and heat transfer as studied by Dayyan et al. [24]. Calculations have been carried out for several values of all the above-mentioned parameters. The present findings show good agreement with previous studies [24], and hence prove their validity. In order to elaborate on this work, some representative outcomes are presented in tabular as well as in plot forms. Tables 1 and 2, respectively, present results for temperature function θ(η) with isothermal wall temperature and isoflux boundary conditions. The temperature function θ(η), decreases for micropolar fluid in comparison with that of Newtonian fluids. It is because that increase in the vortex viscosity of micropolar fluids provides a resistance to the thermal distribution. Table 3, demonstrates the impact of Reynolds number R on − f (0), −θ (0) and −N (0). Table 4 represents the iterative computations of temperature for three the grid sizes and its extrapolated values. Table 5 is presented to exhibit the efficiency of the iterative procedure, it enumerates the number of iterations and ω opt for the grid sizes h, h/2, h/4 when the parameter λ is varied. The increase in R (1 ≤ R ≤ 5) marked an increase in the magnitude of these three physical quantities. These results show correspondence with previous studies by Dayyan et al. [24]. Moreover, the skin friction coefficient − f (0) is lesser in magnitude for micropolar fluids than that of Newtonian fluids; this outcome is supported by the experimental results of Fabula and Hoyt [45], but the heat transfer rate −θ (0) at the surface is greater in magnitude for micropolar fluids than that of Newtonian fluids.   Table 3. The results of skin friction f (0), heat transfer rate θ(0), couple stress N (0) when λ = 0, C 1 = 0.5, Pr = 1 and s = 0 (isothermal).

Micropolar Fluids
Micropolar Fluids  Figure 2 illustrates the reduction in flow speed f (η) with increments in R. Figure 3 shows that the flow speed f (η) is increased when injection parameter λ (λ < 0) increases, this result corresponds to physical nature of the injection/suction parameter. The increasing values of the micropolar parameter C 1 , marked as significant in f (η) and increment in boundary layer thickness as well as shown Figure 4. This result attracts our special attention, as it reveals the clear difference between the dynamics of Newtonian fluids and micropolar fluids. Here, the curve for C 1 = 0 corresponds to Newtonian fluid flow.     Figure 3. Graph of f (isothermal) with variation of λ. Figures 5 and 6 are presented, respectively, to display the non-dimensional microrotation N(η) under the effects of parameters R and λ. It is revealed that N(η) increases near the boundary and then swings down away from the boundary when R and λ (λ > 0) are increased. Figure 7 shows that the increment in the values of micropolar material parameter C 1 shows a notable rise in the curve of N(η). It is reasonable from the physical point of view, that the increment in vortex viscosity causes an increase in the microrotation of the microstructures suspended in the flow field.       and a reduction in thermal distribution as revealed in Figure 8. Furthermore, Figure 9 shows that the   and a reduction in thermal distribution as revealed in Figure 8. Furthermore, Figure 9 shows that the   and a reduction in thermal distribution as revealed in Figure 8. Furthermore, Figure 9 shows that the  Figures 8-11 demonstrate the decreasing pattern of temperature function θ(η), with increasing values of the parameters Pr, λ, C 1 and n , respectively (for isothermal wall temperature boundary conditions). The increase in the value of Pr, means a decrease in thermal conductivity and a reduction in thermal distribution as revealed in Figure 8. Furthermore, Figure 9 shows that the injection λ (λ < 0) raises the temperature function θ(η), but suction λ (λ > 0) causes a reduction in θ(η). The increase in heat index parameter s (0 ≤ s ≤ 10) made a marked reduction in temperature distribution and, hence, a decrease in the thermal boundary layer is observed as depicted in Figure 10. Furthermore, Figure 11 shows that the temperature function θ(η) becomes less for micropolar fluids than that of Newtonian fluids. This is because of the pressure of spinning micro structures in micropolar fluids provide resistance to heat distribution in the bulk fluid.            The results of θ(η), for isoflux boundary conditions have been computed and presented, respectively, in Figures 12-15 under the influence of Pr, R, λ and C 1 . The increment in Pr, R, indicates a reduction in θ(η), respectively, in Figures 12 and 13. This is because of Pr , which is reciprocal to thermal conductivity, and R is reciprocal to the thermal porosity of the medium. Figure 14 illustrates that an increment in injection λ (λ < 0) causes a decrease in θ(η), but an increment in suction at the wall λ (λ > 0) raises the curve of θ(η). The heat function θ(η) is greater for Newtonian fluids than that of micropolar fluids as depicted in Figure 15. indicates a reduction in ) (η θ , respectively, in Figures 12 and 13. This is because of Pr , which is reciprocal to thermal conductivity, and R is reciprocal to the thermal porosity of the medium.    indicates a reduction in ) (η θ , respectively, in Figures 12 and 13. This is because of Pr , which is reciprocal to thermal conductivity, and R is reciprocal to the thermal porosity of the medium.

Conclusions
The free convection boundary layer flow of micropolar fluids through porous medium over a permeable stretching sheet with two thermal boundary conditions which are either power law heat flux or wall temperature is analyzed numerically. The theoretical inclusion of spinning of a suspended microstructure in the bulk flow requires the coupling of angular momentum with linear momentum. Similarly, both the first and second laws of thermodynamics have been applied in heat transfer analysis. A reliable, efficient numerical scheme based on finite difference approximation has been applied with improved accuracy to be achieved through the extrapolation technique. The novelty of the present results included a reduction in flow speed, temperature distribution, skin friction coefficient, and an increase in the heat transfer rate. Moreover, the output in this work presents the results for microrotation and couple stress. The main outcomes of this research work are as follows: The velocity in the boundary layer is retorted against enhanced values of the Reynold's number R , as well as that of suction 0 < λ . In short, the present findings embrace microrotation and micropolar parameters to reveal the

Conclusions
The free convection boundary layer flow of micropolar fluids through porous medium over a permeable stretching sheet with two thermal boundary conditions which are either power law heat flux or wall temperature is analyzed numerically. The theoretical inclusion of spinning of a suspended microstructure in the bulk flow requires the coupling of angular momentum with linear momentum. Similarly, both the first and second laws of thermodynamics have been applied in heat transfer analysis. A reliable, efficient numerical scheme based on finite difference approximation has been applied with improved accuracy to be achieved through the extrapolation technique. The novelty of the present results included a reduction in flow speed, temperature distribution, skin friction coefficient, and an increase in the heat transfer rate. Moreover, the output in this work presents the results for microrotation and couple stress. The main outcomes of this research work are as follows: The velocity in the boundary layer is retorted against enhanced values of the Reynold's number R , as well as that of suction λ < 0 .
The velocity is incremented and the trend of flow becomes faster, increasing the strength of micropolar parameter C 1 or injection (λ > 0). The microrotation is incremented with higher values of C 1 . The microrotation exhibits an oscillatory pattern as it increases near the surface and then decreases away from the surface with increasing values of R and λ.
The temperature function θ(η) and thermal boundary layer thickness decrease against Pr, R, λ, s and C 1 for both cases; isothermal, isoflux, and wall temperature conditions.
Skin friction coefficient − f (0) and local Nusselt number −θ (0) are smaller in value for micropolar fluids than for Newtonian fluids.
Couple stress −N (0) increases with an increase in C 1 .
In short, the present findings embrace microrotation and micropolar parameters to reveal the real nature of fluid motion and thermal distribution. The varying values of micropolar parameter C 1 indicate a notable impact on the physical nature of flow and temperature transportation. In this sense, the present work is more practicable than the previous studies for Newtonian fluids. The numerical scheme used herein is simple and efficient but it may not work for heat and mass transfer across any complex geometry. The numerical solution for the basic partial differential equation may provide more versatile results than the transformed model.  Acknowledgments: The authors are thankful to NTU, Singapore, Mathematics Department, GC University Lahore and UOH, KSA for providing research environment.

Conflicts of Interest:
The authors declare no conflict of interest.
Availability of Data and Material: All data is available and no need of any permission etc.