**Numerical Analysis of Thermal Radiative Maxwell Nanofluid Flow Over-Stretching Porous Rotating Disk**

**Shuang-Shuang Zhou <sup>1</sup> , Muhammad Bilal 2,\* , Muhammad Altaf Khan <sup>3</sup> and Taseer Muhammad 4,5**


**Abstract:** The fluid flow over a rotating disk is critically important due to its application in a broad spectrum of industries and engineering and scientific fields. In this article, the traditional swirling flow of Von Karman is optimized for Maxwell fluid over a porous spinning disc with a consistent suction/injection effect. Buongiorno's model, which incorporates the effect of both thermophoresis and Brownian motion, describes the Maxwell nanofluid nature. The dimensionless system of ordinary differential equations (ODEs) has been diminished from the system of modeled equations through a proper transformation framework. Which is numerically computed with the bvp4c method and for validity purposes, the results are compared with the RK4 technique. The effect of mathematical abstractions on velocity, energy, concentration, and magnetic power is sketched and debated. It is perceived that the mass transmission significantly rises with the thermophoresis parameter, while the velocities in angular and radial directions are reducing with enlarging of the viscosity parameter. Further, the influences of thermal radiation Rd and Brownian motion parameters are particularly more valuable to enhance fluid temperature. The fluid velocity is reduced by the action of suction effects. The suction effect grips the fluid particles towards the pores of the disk, which causes the momentum boundary layer reduction.

**Keywords:** bvp4c; RK4 technique; brownian motion; porous rotating disk; maxwell nanofluid; thermally radiative fluid; von karman transformation

**1. Introduction**

The researchers have been interested in Maxwell nanofluid flow over a porous spinning disc because of its many uses in engineering and innovation. Non-Newtonian fluids are important in a variety of manufactured liquids, including plastics, polymers, pulps, toothpaste and fossil fluids. To simulate the analysis of these liquids, a variety of models have been suggested. Shear stress and shear rate are linked in non-Newtonian liquids because of their nonlinear existence. The momentum equation in these fluids involves dynamic nonlinear terms, making it difficult to solve. A variety of mathematical models exist in the literature to simulate the performance of these fluids.

In the present era, the role of nanotechnology to fulfill the increasing demand for energy and face energy challenges is remarkable. The usage of nanoparticles in ordinary base fluid (water, kerosene oil, etc.) effectively enhances the heat transfer and improves their thermal properties. The applications of nanoparticles in certain fields of engineering and industry are in the cooling systems of electronic devices and in cancer therapy, heat

**Citation:** Zhou, S.-S.; Bilal, M.; Khan, M.A.; Muhammad, T. Numerical Analysis of Thermal Radiative Maxwell Nanofluid Flow Over-Stretching Porous Rotating Disk. *Micromachines* **2021**, *12*, 540. https:// doi.org/10.3390/mi12050540

Academic Editors: Lanju Mei and Shizhi Qian

Received: 21 March 2021 Accepted: 3 May 2021 Published: 10 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

exchangers, transformer cooling, nuclear reactors, and space cooling systems. Due to the ability of oil wetting and dispersing, they are also used for cleaning purposes, in power generation, microfabrication, hyperthermia, and metallurgical purposes.

#### **2. Literature Review**

The rotating disk phenomena are widely used in centrifugal filtration, turbomachines, the braking system of vehicles, jet motors, sewing machines, turbine systems, heat exchangers and computer disk drives, etc. Von Karman [1] for the first time introduced similarity transformation. To solve Navier-Stokes equations, he studied fluid flow over an infinite rotating disk. Cochran [2] employed Von Karman's similarity transformation to incompressible fluid over a rotating disk and examined the asymptotic solution. Wagner [3] investigated the mechanism of heat transfer over the rotating disk, by considering Von Karman's velocity distribution, and analyzed convection in the non-turbulent flow. Turkyimazogl [4] looked at fluid movement over a spinning disc that was stretching under the influence of a static electric field. Liang et al. [5] reported a comparative study between semi-analytical model and experimental data to yield the best settlement. Millsaps and Pohlhausen [6] used Von Karman's similarity approach and analyzed heat distribution with the consequences of entropy generation over revolving disk. The three-dimensional (3D) magnetohydrodynamics (MHD) stagnation flow of ferrofluid, the numerical solution was revealed by Mustafa et al. [7]. Mustafa et al. [8], by taking MHD nanofluid over rotating surface with the effects of partial slips, observed that the boundary layer thickness and momentum transport are reduced due to slip effects. Rashidi et al. [9] have used a spinning disc to perform a viscous dissipation review for MHD nanofluid.

The people of the modern world are facing many challenges due to the increasing demand for energy by the latest technologies. Firstly Choi [10] presented the nanofluids terminology. The Brownian motion and thermophoresis mechanisms bring about a significant role in improving the thermal properties of base fluid presented by Buongiorno [11]. Turkyilmazolglu [12] analytically studied the energy and momentum equations of nanofluid flow, to deduce heat and flow transport. Pourmehran [13] considered Cu and Al2O<sup>3</sup> nanoparticles to study heat and flow transfer in the microchannel. Hatami et al. [14] reported the heat transfer in nanofluid with the phenomena of natural convection. The Oldroyd-B fluid with nanoparticles over stretching sheet surface was reported by Nadeem et al. [15]. Aziz and Afify [16] have used the technique of the Lie group, to study non-Newtonian nanofluids. Yang et al. [17] studied the convective heat with Buongiorno Model's for nanofluid in the concentric annulus.

The study of a Newtonian fluid, due to its wide applicability in different fields of science and engineering, attracted the attention of scientists and researchers during the last few decades. Its major role is in geophysics, polymer solution, paper production, cosmetic processes, exotic lubricants, paints, suspensions, colloidal solutions, nuclear and chemical industries, pharmaceuticals, oil reservoirs, bioengineering, etc. [18]. Xiao et al. [19] attempted a fractal model for the capillary flow through a torturous capillary with the non-smooth surface in porous media. Attia [20] evaluated the Reiner-Rivlin numerical simulations for thermal convection over a porous spinning disc qualitatively. Griffiths [21] tested the Newtonian fluid Carreau viscosity model and high shear stresses on spinning discs. The numerical analysis of Reiner-Rivlin fluid flow for heat transfer and slip flow over a spinning disc is treated by Mustafa and Tabassum [22]. The micropolar fluid for thermophoretic diffusion generated by the rotation of the disk was examined by Doh and Muthtamilselvan [23].

Darcy's law is a mathematical equation that explains how fluid flows through a porous medium. Henry Darcy developed the law based on the effects of studies on the flow of water across sand beds, laying the groundwork for hydrogeology, a branch of earth sciences [24]. Fourier's law in heat conduction, Ohm's law in electrical networks, and Fick's law in diffusion theory are all examples of this law. Morris Muskat [25] improved Darcy's equation for a single-phase flow by incorporating viscosity into Darcy's single (fluid) phase equation. It is easy to see that viscous fluids have a harder time passing through a porous medium than less viscous fluids. Rasool et al. [26–28] numerically simulated the Darcy-Forchheimer effect on MHD nanoliquid flow between stretching non-linear sheets. Rasool et al. [29] scrutinized the consequences of thermal radiation, chemical reaction and Dufour-Soret on incompressible steady Darcy-Forchheimer flow of nanoliquid. Shafiq et al. [30] studied nanofluid flow under the influence of convective boundary conditions and thermal slip over a spinning frame. They found that the axial and transverse velocity fields all drop significantly due to the Forchheimer number's strong retardation. Skin friction is intensified by the Forchheimer number and porosity ratios, while skin friction is diminished by all slip parameters.

Viscoelastic fluids are a subclass of Newtonian fluid having memory effects. The intensity of energy discharged by these fluids is mainly accountable for recovery after the stress is removed. The Maxwell flow regime is the most basic viscoelastic fluid model, expressing memory effects by fluid relaxation time [31]. The attitude of the current model is very close to that of other geomaterials and polymers models. The aim of the present work is to provide a mathematical model for unsteady boundary layer flow of non-Newtonian Maxwell nanofluid with the heat transmission over a porous spinning disc. The present work has many industrial and engineering applications, which increases its worth. Using a resemblance method, the system of ODEs is limited to a structure of PDEs. A boundary value solver (bvp4c) technique is used to draw a numerical solution to the problem while RK4 method has been applied for validity.

#### **3. Formulation of the Problem**

Consider an unsteady hybrid nanoliquid flow over a stretching porous spinning disc. The magnetic force *B*<sup>0</sup> is introduced to the disc vertically. The disc rotates and stretches at different speeds (*u*, *v*) = (*cr*, *c*Ω), where *c* and *w* are the spinning and extending rates, respectively. The disk temperature is represented by *τw*. The formulation of the problems is conducted in (*r*, *ϕ*, *z*) cylindrical coordinates, where *u*, *v*, *w* is velocity component increasing in (*r*, *ϕ*, *z*) direction. At *z*− the axis, the motion of the disk is assumed to be axisymmetric. The thermal radiation is significant in modeling the energy equation. The viscosity of a fluid is taken to be temperature-dependent *µ*(*τ*) = *µ*0*e* −*ζ*(*τ*−*τ*0) . The concentration and temperature are represented by (*Cw*, *τw*) and (*C*∞, *τ*∞) represent the concentration and temperature above the disk surface.

#### *3.1. Governing Equations*

Under the presuppositions stated above, the flow equations are as observes [31]:

$$
\nabla.\mathbf{V} = \mathbf{0},
\tag{1}
$$

$$
\rho\_f(\mathbf{V}.\nabla)\mathbf{V} = \nabla p + \nabla.\mathbf{S} + j \times \mathbf{B},\tag{2}
$$

$$(V.\nabla)\tau = a\nabla^2\tau + \tau^\* \left( D\_B \nabla \mathbb{C}.\nabla \tau + \frac{D\tau}{\tau\_{\infty}} \nabla \tau. \nabla \tau \right) + \nabla.q\_{rad,} \tag{3}$$

$$(V.\nabla)\mathcal{C} = D\_B\nabla^2\mathcal{C} + \frac{D\tau}{\tau\_\infty}\nabla^2\tau,\tag{4}$$

$$
\nabla.B = 0,\tag{5}
$$

$$
\rho \frac{\partial B}{\partial t} = \rho \nabla + (V \times \mathcal{B}) + \frac{\rho}{\sigma \mu\_2} \nabla^2 V. \tag{6}
$$

where *DB*, *DT*, *ρ<sup>f</sup>* , *V* and *α* are the coefficients of Brownian motion, thermophoretic diffusion, fluid density, velocity, and thermal and diffusivity respectively. Equations (1)–(6) are simplified because of the boundary layer approximation concept:

$$
\frac{\partial u}{\partial r} + \frac{u}{r} + \frac{\partial w}{\partial z} = 0,\tag{7}
$$

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} - \frac{v^2}{r} = \left[ \begin{array}{c} \frac{1}{\rho} \left( \mu(r) \frac{\partial u}{\partial z} \right) - \lambda\_1 (u^2 \frac{\partial^2 u}{\partial r^2} + w \frac{\partial^2 u}{\partial z^2} + 2uw \frac{\partial^2 u}{\partial r \partial z} - \frac{2uv}{r} \right) \\\ \frac{\partial v}{\partial r} - \frac{2vw}{r} \frac{\partial v}{\partial z} + \frac{uv^2}{r^2} + \frac{v^2}{r^2} \frac{\partial u}{\partial r} \end{array} \right], \tag{8}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} - \frac{uv}{r} = \begin{bmatrix} \frac{1}{\rho} \left( \mu(\tau) \frac{\partial v}{\partial z} \right) - \lambda\_1 (u^2 \frac{\partial^2 v}{\partial r^2} + w^2 \frac{\partial^2 v}{\partial z^2} + 2uw \frac{\partial^2 v}{\partial r \partial z} + \frac{2uv}{r} \right) \\\ \frac{\partial u}{\partial r} + \frac{2vw}{r} \frac{\partial u}{\partial z} + \frac{2u^2 v}{r^2} - \frac{v^2}{r^2} \frac{\partial v}{\partial r} \end{bmatrix}, \tag{9}$$

$$\frac{\partial \tau}{\partial t} + \mu \frac{\partial \tau}{\partial r} + w \frac{\partial \tau}{\partial z} = \frac{k}{\rho\_{c\_p}} \left( \frac{\partial^2 \tau}{\partial z^2} \right) + \tau^\* (D\_B \frac{\partial \tau}{\partial z} \frac{\partial C}{\partial z} + \frac{D\_\tau}{\tau\_\infty} \left( \frac{\partial \tau}{\partial z} \right)^2) - \frac{1}{\rho\_{c\_p}} \frac{\partial q\_r}{\partial z}, \tag{10}$$

$$\frac{\partial \mathcal{C}}{\partial t} + \mu \frac{\partial \mathcal{C}}{\partial r} + w \frac{\partial \mathcal{C}}{\partial z} = D\_B \left( \frac{\partial^2 \mathcal{C}}{\partial z^2} \right) + \frac{D\_\tau}{\tau\_\infty} \left( \frac{\partial^2 \tau}{\partial z^2} \right), \tag{11}$$

$$\frac{\partial B\_{\rm r}}{\partial t} = \left[ -w \frac{\partial B\_{\rm r}}{\partial z} - B\_{\rm r} \frac{\partial w}{\partial z} + u \frac{\partial B\_{\rm z}}{\partial z} + B\_{\rm z} \frac{\partial u}{\partial z} + \frac{1}{\sigma \mu\_2} \left( \frac{\partial^2 B\_{\rm r}}{\partial r^2} + \frac{\partial^2 B\_{\rm r}}{\partial z^2} + \frac{1}{r} \frac{\partial B\_{\rm r}}{\partial r} - \frac{B\_{\rm r}}{r^2} \right) \right], \tag{12}$$

$$\frac{\partial B\_{\overline{z}}}{\partial t} = \left[ w \frac{\partial B\_{r}}{\partial r} + B\_{r} \frac{\partial w}{\partial r} + \frac{1}{r} w B\_{r} - u \frac{\partial B\_{\overline{z}}}{\partial r} - B\_{\overline{z}} \frac{\partial u}{\partial r} - \frac{1}{r} u B\_{\overline{z}} + \frac{1}{\sigma \mu\_{2}} \left( \frac{\partial^{2} B\_{\overline{z}}}{\partial r^{2}} + \frac{\partial^{2} B\_{\overline{z}}}{\partial z^{2}} + \frac{1}{r} \frac{\partial B\_{\overline{z}}}{\partial r} \right) \right]. \tag{13}$$

The temperature difference within a flow is assumed to be small, therefore higherorder in Taylor series and ignored at *τ*∞, by using Rosseland approximation, the simple form of radiation heat flux is as follow [31]:

$$q\_r = \frac{-4\sigma^\*\partial\tau^4}{3k^\*\partial z} = -\frac{16\sigma^\*\tau^3}{3k^\*}\frac{\partial\tau}{\partial z'}\tag{14}$$

using Equation (14) in (10), we get

$$\frac{\partial \tau}{\partial t} + u \frac{\partial \tau}{\partial r} + w \frac{\partial \tau}{\partial z} = \frac{k}{\rho\_{\varepsilon\_p}} \left( \frac{\partial^2 \tau}{\partial z^2} \right) + \tau^\* \left( D\_B \frac{\partial \tau}{\partial z} \frac{\partial \mathcal{C}}{\partial z} + \frac{D\_\tau}{\tau\_\infty} (\frac{\partial \tau}{\partial z})^2 \right) - \frac{16 \sigma^\* \tau^3}{3k^\*} \frac{\partial \tau}{\partial z}. \tag{15}$$

The boundary conditions are:

$$\begin{array}{ccccc} \mu = cr, v = \Omega r, \; w = W, \; \mathbb{C} = \mathbb{C}\_{w}, \; T = T\_{w}, \; B\_{r} = 1, B\_{\mathbb{Z}} = 1 & at \; z = 0 \\\ u \to 0, \; v \to 0, \; w \to 0, \; \mathbb{C} \to \mathbb{C}\_{\infty}, \; T \to T\_{\infty}, \; B\_{r} \to 0, B\_{\mathbb{Z}} \to 0 \; as & Z \to \infty. \end{array} \tag{16}$$

## *3.2. Similarity Transformation*

Considering the following transformation, to reduce the system of PDEs to the system of ODEs:

$$\begin{aligned} \boldsymbol{u} &= \frac{\boldsymbol{c}\boldsymbol{r}}{1 - \boldsymbol{a}t} \boldsymbol{F}(\eta), \boldsymbol{\nu} = \frac{\Omega \boldsymbol{r}}{1 - \boldsymbol{a}t} \boldsymbol{G}(\eta), \; \boldsymbol{w} = \frac{\sqrt{\boldsymbol{c}\boldsymbol{v}}}{1 - \boldsymbol{a}t} \boldsymbol{H}(\eta), \boldsymbol{\eta} = \sqrt{\frac{\boldsymbol{c}}{\boldsymbol{v}}} \boldsymbol{z}, \; \mathcal{B}\_{\boldsymbol{r}} &= \frac{\boldsymbol{c}\boldsymbol{M}\_{0}}{1 - \boldsymbol{a}t} \boldsymbol{M}(\eta), \\ \boldsymbol{B}\_{\boldsymbol{z}} &= \frac{-\boldsymbol{M}(2\boldsymbol{v}\_{f}\boldsymbol{\epsilon})^{\frac{1}{2}}}{1 - \boldsymbol{a}t} \boldsymbol{N}(\eta), \; \; \boldsymbol{T} = (T\_{\infty}) + \boldsymbol{\Theta}(\eta)(T\_{W} - T\_{\infty}) \boldsymbol{C} = (\mathbb{C}\_{\infty}) + \boldsymbol{\phi}(\eta)(\mathbb{C}\_{W} - \mathbb{C}\_{\infty}). \end{aligned} \tag{17}$$

The following system of ODEs is obtained by using Equation (17) in Equations (7)–(13) and (15) and (16):

$$F^{r} = F^{\ell} \Theta^{\ell} + \frac{\varepsilon^{\theta \Theta}}{\delta} \left\{ S \left( \frac{F^{\prime} \eta}{2} + F \right) + F^{2} + HF^{\prime} - G^{2} \right\} + \beta\_{1} \frac{\varepsilon^{\theta \Theta}}{\delta} \left\{ HF^{\prime} + 2FF^{\prime}H - 2HG^{\prime} \right\} - M \frac{\varepsilon^{\theta \Theta}}{\delta} \left( \mathcal{F} + \beta\_{1} HF^{\prime} \right), \tag{18}$$

$$G'' = \frac{\delta G' \Theta' + \epsilon^{\delta \Theta} \left\{ S \left( \frac{G' - \eta}{2} + G \right) + 2FG + HG' \right\} + \beta\_1 \epsilon^{\delta \Theta} \left\{ 2FHG' + 2FHG \right\} + Me^{\delta \Theta} (G + \beta\_1 HG')}{1 + \beta\_1 H^2 \epsilon^{\delta \Theta}},\tag{19}$$

$$\Theta'' = \frac{-4Rd\Theta^2\Theta^2(\Theta - 1)^3 - 6\Theta'^2\Theta(Q\_w - 1)^2 - 3\Theta'^2(Q\_w - 1) - \Pr\left(-\frac{5}{2}(\Theta'\eta) - R\Theta' + Nb\Theta\phi' + Nb\Theta^3\right)}{\left(1 + \frac{4}{3}\right)\mathcal{R}d + \frac{4}{3}\mathcal{R}d\Theta^3(Q\_w - 1)^3 + 3\Theta^2(Q\_w - 1)^2 + 3\Theta(\Theta - 1)}\,,\tag{20}$$

$$
\Phi'' = \mathcal{Sc}\left(A\phi'\eta + H\phi'\right) - \frac{\mathcal{N}t}{\mathcal{N}b}\Theta'',\tag{21}
$$

$$M''' = Bt \left[ -HM'' + M'H' + FN' + NF' + S \left( \frac{M''\eta}{2} + M' \right) \right],\tag{22}$$

$$M'' = -Bt\left[2HM' + 2NF - \frac{S}{2}(N\eta + N)\right].\tag{23}$$

the transforms conditions are:

$$\begin{array}{l} F(0) = 1, \ G(0) = \omega, \ H = W\_{\sf s}, \ \Theta(0) = 1, \ f(0) = 1, \ M'(0) = 0, \ N(0) = 1, \\\ F(\infty) = 0, \ G(\infty) = 0, \ H(0) = 0, \ \Theta(\infty) = 0, \ f(\infty) = 0, \ M(\infty) = 0, \ N(\infty) = 0. \end{array} \tag{24}$$

where *Pr* = *ν*/*α* is the Prandtl number, *Sc* = *v*/*D<sup>B</sup>* is the Schmidt number, *β*<sup>1</sup> = *λ*1*c* the Deborah number, while suction/injection parameter, Brownian motion, variable viscosity, thermophoresis parameter, magnetic parameter, temperature ratio, and thermal radiation are defined as [31]:

$$\text{Ws} = \frac{w}{\sqrt{\text{c}\mathcal{V}}}, \text{N} \\ b = \frac{\text{\textdegree } ^\text{s}D\_{\text{B}}(\text{C}\_{\text{w}} - \text{C}\_{\text{ox}})}{v}, \delta = \text{\textdegree } (\text{\textdegree } \pi\_{\text{w}} - \text{\textdegree } \pi\_{\text{c}}), \text{N} \\ t = \frac{\text{\textdegree } ^\text{s}D\_{\text{\textdegree }}(\pi\_{\text{w}} - \pi\_{\text{c}})}{\pi\_{\text{c}}v}, \text{M} = \frac{\sigma B\_{0}^{2}}{c\_{p}}, \text{R} \\ d = \frac{4 \sigma^{\*} T\_{\infty}^{3}}{k k^{\*}}. \tag{25}$$

The skin friction *C fx*, local Sherwood *Sh<sup>r</sup>* and Nusselt number *Nu<sup>r</sup>* are mathematically can be written as [32,33]:

$$\mathcal{C}f = \frac{\sqrt{\pi\_{zr}^{2} +} \pi\_{z\rho}^{2}}{\rho \left(\Omega r\right)^{2}},\tag{26}$$

$$Sh\_{\prime} = -\frac{r}{\left(\mathcal{C}\_{w} - \mathcal{C}\_{\infty}\right)} \left(\frac{\partial \tau}{\partial \mathcal{C}}\right)\Big|\_{z=0\prime} \tag{27}$$

$$Nu\_r = \frac{r}{\left(\tau\_w - \tau\_\infty\right)} \left[1 + \frac{16\sigma^\*\tau^3}{3kk^\*}\right] \left(\frac{\partial \tau}{\partial z}\right)\Big|\_{z=0}.\tag{28}$$

The skin friction, Sherwood and Nusselt numbers have a non-dimensional structure as:

$$\mathbb{C}f\text{Re}\_{\mathbf{r}}^{1/2} = \sqrt{f''^2(0) + \mathbf{g}'^2(0)},\tag{29}$$

$$Re^{-}\frac{1}{2}Sh\_{\!\!\!r} = -\phi(0)\_{\!\!\!\!r} \tag{30}$$

$$Re^{-}\frac{1}{2}Nu\_{l} = -\left(1 + \frac{4}{3}Rd\{1 + (\Theta\_{\overline{w}} - 1)\Theta(0)\}^3\right)\Theta'(0).\tag{31}$$

Here *Re* = *ru ν* is the local Reynold number.

#### **4. Solution Procedures**

The higher-order model equation is brought down to first order by choosing variables:

$$\begin{aligned} \chi\_1 &= H, \; \chi\_2 = F, \; \chi\_3 = F\prime, \; \chi\_4 = G\prime, \; \chi\_5 = G\prime, \; \chi\_6 = \theta\prime, \; \chi\_7 = \theta\prime, \\\chi\_8 &= \phi\prime, \; \chi\_9 = \phi\prime, \; \chi\_{10} = M\prime, \; \chi\_{11} = M\prime, \; \chi\_{13} = N\prime, \; \chi\_{14} = N\prime. \end{aligned} \tag{32}$$

*χ* 0 <sup>1</sup> = −2*χ*2, *χ* 0 <sup>1</sup> = *χ*3, *χ* 0 <sup>3</sup> = *<sup>χ</sup>*7*χ*<sup>3</sup> + *<sup>e</sup> δθ δ S*( *ηχ*3 <sup>2</sup> + *χ*2) + *χ* 2 <sup>2</sup> + *χ*1*χ* 2 <sup>2</sup> − *χ* 2 4 + *β*<sup>1</sup> *e δθ δ* {*χ*1*χ*<sup>3</sup> + 2*χ*1*χ*2*χ*<sup>3</sup> − 2*χ*1*χ*4*χ*5} *M*<sup>0</sup> *e δθ δ* (*χ*<sup>2</sup> + *β*1*χ*1*χ*3), *χ* 0 <sup>4</sup> = *χ*5, *χ* 0 <sup>5</sup> = *δχ*5*χ*7+*e δθ*n *S*( *χ*5−*η* 2 )+2*χ*2*χ*4+*χ*1*χ*<sup>5</sup> o +*β*1{2*χ*1*χ*2*χ*5−2*χ*1*χ*4*χ*2}+*M*0*e δθ* (*χ*4+*β*1*χ*1*χ*5) 1+*β*<sup>1</sup> *e δθχ* 2 1 , *χ* 0 <sup>6</sup> = *χ*7, *χ* 0 <sup>7</sup> = 4 3 *Rdχ* 2 7 (*θw*−1) n 3*χ* 2 6 (*θw*−1) <sup>2</sup>+6*χ*6(*θw*−1)−<sup>3</sup> o <sup>−</sup>Pr{*Nbχ*7*χ*9+*Nt<sup>χ</sup>* 2 <sup>7</sup>−*Aχ*7*χ*10−*χ*1*χ*7} <sup>1</sup><sup>−</sup> <sup>4</sup> 3 *Rd*<sup>−</sup> <sup>4</sup> 3 *Rdχ*6(*θw*−1){(*θw*−1)*<sup>χ</sup>* 2 <sup>6</sup>−3*χ*6(*θw*−1)+3} , *χ* 0 <sup>8</sup> = *χ*9,

$$\chi\_{9}^{\prime} = \frac{\text{Nt}}{\text{Nb}} \left( \frac{\frac{4}{3} \text{R} \chi\_{7}^{2} (\theta\_{\text{w}} - 1) \left\{ 3 \chi\_{6}^{2} (\theta\_{\text{w}} - 1)^{2} + 6 \chi\_{6} (\theta\_{\text{w}} - 1) - 3 \right\} - \text{Pr} \left\{ \text{Nb} \chi\_{7} \chi\_{8} + \text{Nt} \chi\_{7}^{2} - \text{A} \chi \chi\_{10} - \text{X} \chi \chi\_{1} \right\}}{1 - \frac{4}{3} \text{R} \text{d} - \frac{4}{3} \text{R} \text{d} \chi\_{6} (\theta\_{\text{w}} - 1) \left\{ (\theta\_{\text{w}} - 1) \chi\_{6}^{2} - 3 \chi\_{6} (\theta\_{\text{w}} - 1) + 3 \right\}} \right)$$
  $\text{Sc} (A \chi\_{9} \chi\_{10} + \chi\_{9} \chi\_{9})$   $\chi\_{1} = \text{Sc} (A \chi\_{9} \chi\_{10} + \chi\_{9} \chi\_{9})$   $\chi\_{1} = \text{Re} \left[ -2 \chi\_{11} \chi\_{2} - \chi\_{1} \chi\_{11} + \chi\_{2} \chi\_{14} + \chi\_{13} \chi\_{2} + \text{S} \left( \frac{\eta\_{1} \chi\_{12}}{2} + \chi\_{11} \right) \right]$   $\chi\_{1}^{\prime} = \chi\_{11} \chi\_{1}^{\prime} - \text{Re} \left[ 2 \chi\_{11} \chi\_{1} + 2 \chi\_{13} \chi\_{2} + \frac{\xi}{2} (\eta \chi\_{14} + \chi\_{13}) \right]$ 

The boundary conditions are:

$$\begin{aligned} \chi\_1(0) &= 1, \chi\_2(0) = 0, \chi\_4(0) = 1, \chi\_6(0) = 1, \chi\_8(0) = 1, \chi\_{11}(0) = 0, \chi\_{13}(0) = 1, \\\chi\_2(\infty) &= 0, \chi\_4(\infty) = 0, \chi\_6(\infty) = 0, \chi\_8(\infty) = 0, \chi\_{11}(\infty) = 1, \chi\_{13}(\infty) = 0. \end{aligned} \tag{34}$$

#### **5. Results and Discussion**

The discussion section is devoted to understanding better the graphical and physical description. The system of non-linear Equations (18)–(23) along with their boundary conditions, Equation (24), are solved through numerical method bvp4c. The configuration of the problem is described in Figure 1. The velocities, energy profile, concentration distribution *φ*(*η*), magnetic strength in the radial direction *M*(*η*), and azimuthal magnetic strength *N*(*η*) are explored graphically through different physical constraints Figures 2–10. While keeping *Pr* = 6.7, *ω* = 1.0, *δ* = 0.5, *β* = 0.1, *M* = 1.1, *Nb* = 0.5, *Nt* = 0.7, *Q<sup>w</sup>* = 1.2, *Sc* = 2.0, and *Rd* = 0.5. *Micromachines* **2021**, *12*, x FOR PEER REVIEW 7 of 15 **5. Results and Discussion**  The discussion section is devoted to understanding better the graphical and physical description. The system of non-linear Equations (18)–(23) along with their boundary conditions, Equation (24), are solved through numerical method bvp4c. The configura-

Figure 2a–c depicts the behavior velocities profiles against the variation of Deborah number *β*. All three velocity shows decreasing behavior for incremented of *β*. Deborah number is the measure of the content evaluation period to content recreational time, so having optimum stress relaxation or eliminating observation time increases the value of *β*. It reflects the fluid's solid-like reaction. The hydrodynamic boundary layer thins out, and the velocity experiences more resistance. tion of the problem is described in Figure 1. The velocities, energy profile, concentration distribution φ ( ) η , magnetic strength in the radial direction *M* ( ) η , and azimuthal magnetic strength *N* ( ) η are explored graphically through different physical constraints Figures 2–10. While keeping *Pr* = 6.7 , ω = 1.0 , δ = 0.5 , β = 0.1 , *M* = 1.1 , *Nb* = 0.5 , *Nt* = 0.7 , 1.2 *Qw* = , *Sc* = 2.0 , and *Rd* = 0.5 .

Figure 2a–c depicts the behavior velocities profiles against the variation of Deborah

β. Debo-

. All three velocity shows decreasing behavior for incremented of

rah number is the measure of the content evaluation period to content recreational time, so having optimum stress relaxation or eliminating observation time increases the value

. It reflects the fluid's solid-like reaction. The hydrodynamic boundary layer thins

The implications of the suction factor ( 0) *Ws* < on the velocities profile and heat spectrum on the disc surface are depicted in Figure 3a–c. It seems to be that as the suction velocity rises, the velocity decreases. Because the suction velocity draws the fluid parti-

Figures 4a–d are drawn to depict the effects of injection velocity parameter ( 0) *Ws* > on the radial, azimuthal, tangential velocity profile and temperature distribution. It can be seen that the enhancing of injection velocity increases the velocity and fluid tempera-

The radial, azimuthal, tangential velocity profile and temperature distribution variation are illustrated in Figure 5a–d. Figure 5a depicts the dominant behavior of radial

shows that the rotation parameter becomes greater than extending. The radial velocity exceeds the disc stretching velocity too close to the disc stretching surface. It does, however, gradually fade away from the disc. It is concluded that the impact of centrifugal force is limited and dominant in the vicinity of the disk's surface. Besides, angular veloc-

. Physically, fluid particles are moved in a radial direction owing to

near the disk flourished shown in Figure 5b. The axial velocity increasing

ω

. The increasing value of

ω

**Figure 1.** Stretchable porous rotating disk.

out, and the velocity experiences more resistance.

**Figure 1.** Stretchable porous rotating disk.

number

of β

ture.

velocity with

ity *G* ( ) η

ω

β

centrifugal force with the enhancement of parameters

respectively.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 9 of 15

**Figure 2.** The out-turn of Deborah number on the axial, radial and azimuthal velocity profiles, **Figure 2.** The out-turn of Deborah number *β* on the axial, radial and azimuthal velocity profiles, respectively. **Figure 2.** The out-turn of Deborah number on the axial, radial and azimuthal velocity profiles, respectively.

**Figure 3.** The out-turn of the suction parameter ( ൏ 0) on the axial velocity, radial velocity and temperature distribution profiles, respectively. **Figure 3.** The out-turn of the suction parameter ( ൏ 0) on the axial velocity, radial velocity and temperature distribution profiles, respectively. **Figure 3.** The out-turn of the suction parameter (*Ws* < 0) on the axial velocity, radial velocity and temperature distribution profiles, respectively.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 10 of 15

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 10 of 15

**Figure 4.** The out-turn of injection parameter ( 0) on the axial, radial and azimuthal velocity profiles and temperature distribution, respectively. **Figure 4.** The out-turn of injection parameter (*Ws* > 0) on the axial, radial and azimuthal velocity profiles and temperature distribution, respectively. **Figure 4.** The out-turn of injection parameter ( 0) on the axial, radial and azimuthal velocity profiles and temperature distribution, respectively.

**Figure 5.** The out-turn of rotation parameter on the axial, radial and azimuthal velocity and **Figure 5.** The out-turn of rotation parameter on the axial, radial and azimuthal velocity and temperature profiles, respectively. **Figure 5.** The out-turn of rotation parameter *<sup>ω</sup>* on the axial, radial and azimuthal velocity and temperature profiles, respectively.

temperature profiles, respectively.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 11 of 15

**Figure 6.** The out-turn of viscosity parameter on the axial, radial and azimuthal velocity and **Figure 6.** The out-turn of viscosity parameter *δ* on the axial, radial and azimuthal velocity and temperature profiles, respectively. **Figure 6.** The out-turn of viscosity parameter on the axial, radial and azimuthal velocity and temperature profiles, respectively.

temperature profiles, respectively.

profile.

profile.

**Figure 7.** The out-turn of different parameters (, , Θ, ) on the temperature distribution **Figure 7.** The out-turn of different parameters (, , Θ, ) on the temperature distribution **Figure 7.** The out-turn of different parameters (*Pr*, *Rd*, Θ, *Nb* ) on the temperature distribution profile.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 12 of 15

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 12 of 15

**Figure 8.** The out-turn of different parameters (, , , ) on the concentration profile. **Figure 8.** The out-turn of different parameters (*Nt*, *Nb*, *Ws*, *Sc* ) on the concentration profile. **Figure 8.** The out-turn of different parameters (, , , ) on the concentration profile.

**Figure 9.** The out-turn of different parameters (, ) on the axial and radial magnetic strength profiles, respectively. **Figure 9.** The out-turn of different parameters (, ) on the axial and radial magnetic strength profiles, respectively. **Figure 9.** The out-turn of different parameters (*Bt*, *Ws* ) on the axial and radial magnetic strength profiles, respectively.

ω

*Rd* = 0.5 .

**Figure 10.** The out-turn of radiation parameter and Brownian coefficient versus the heat and mass **Figure 10.** The out-turn of radiation parameter and Brownian coefficient versus the heat and mass transmission profiles, respectively.

transmission profiles, respectively. **Table 1.** The numerical outcomes for skin friction <sup>ᇱ</sup> (0).  **Mustafa et al. [7] Ahmed et al. [31] Present Paper**  0 −1.1737 −1.1379 −1.1380 The implications of the suction factor (*Ws* < 0) on the velocities profile and heat spectrum on the disc surface are depicted in Figure 3a–c. It seems to be that as the suction velocity rises, the velocity decreases. Because the suction velocity draws the fluid particles towards the pores in the disk, which causes the momentum boundary layer reduction. The enhancement of suction velocity also decreases the fluid temperature.

> 1 −0.9483 −0.9485 −0.9487 2 −0.3262 −0.3264 −0.3266 5 3.1937 3.11937 3.11939 Figure 4a–d are drawn to depict the effects of injection velocity parameter (*Ws* > 0) on the radial, azimuthal, tangential velocity profile and temperature distribution. It can be seen that the enhancing of injection velocity increases the velocity and fluid temperature.

10 12.7209 12.7209 12.7811 20 40.9057 40.9057 40.9058 **Table 2.** Comparison between Runge Kutta order four and bvp4c method. **RK4 bvp4c Absolute Error**  1.0 1.000000 1.000000 8.146310 × 10ିଵଷ 1.2 1.199831 1.199831 3.387821 × 10ିଽ 1.4 0.988459 0.988459 2.845561 × 10ିଽ 1.6 0.879189 0.879189 2.813281 × 10ିଽ 1.8 0.539393 0.539393 3.287961 × 10ିଽ The radial, azimuthal, tangential velocity profile and temperature distribution variation are illustrated in Figure 5a–d. Figure 5a depicts the dominant behavior of radial velocity with *ω*. Physically, fluid particles are moved in a radial direction owing to centrifugal force with the enhancement of parameters *ω*. The increasing value of *ω* shows that the rotation parameter becomes greater than extending. The radial velocity exceeds the disc stretching velocity too close to the disc stretching surface. It does, however, gradually fade away from the disc. It is concluded that the impact of centrifugal force is limited and dominant in the vicinity of the disk's surface. Besides, angular velocity *G*(*η*) near the disk flourished shown in Figure 5b. The axial velocity increasing against the strengthening of *ω* is illustrated in Figure 5c. Figure 5d demonstrates the decreases of temperature *θ*(*η*) with *ω*. Physically, a faster spinning disc reduces the width of the thermal boundary, which play a significant role in the cooling of the system.

**Table 3.** The comparison of *RK* 4 and Bvp4c for Sherwood *<sup>r</sup> Sh* and Nusselt number *Nur* , while keeping ω =1.0 , δ = 0.5 , β = 0.1 , *M* = 1.1, *Nb* = 0.5 , *Nt* = 0.7 , 1.2 *Qw* = and Figure 6a–d are rough sketches that show the effects of the viscosity factor on the velocity and temperature spectrum on the disk's surface. The viscosity parameter causes radial velocity, tangential velocity, and temperature profile to increase, while azimuthal velocity reduces.

*Shr**Nur* **Pr Bvp4c** *RK* **4** *Sc* **Bvp4c** *RK* **4** 3.0 0.7718284 0.7718285 1.0 1.457986 1.457995 4.0 0.6999885 0.6999885 1.5 1.531211 1.531220 5.0 0.6290089 0.6290088 2.0 1.596949 1.596949 6.0 0.5632372 0.5632370 2.5 1.685639 1.685639 **6. Conclusions**  In the present mathematical model, the Maxwell nanoliquid flow over a porous Figure 7a,b indicates how the heat transfer changes when the Prandtl number *Pr* and the thermal radiation factor change. The thermal dispersion of a strong Prandtl fluid is low, while the thermal diffusivity of a low Prandtl fluid is high. Figure 7b is drawn to explore the influence of radiation parameters *Rd* on the thermal mechanism of fluid. With a higher value of the radiation parameter, an increasing trend in heat is analyzed. Because, with enhancing of thermal radiation the fluid absorbs more heat, and as a result increment occurs in boundary layer thickness and fluid temperature. Figure 7c is sketched to illustrate temperature *θ*(*η*) variation versus temperature ratio *Qw*. While Figure 7d indicates the consequence of molecular diffusion on the thermal performance of a Maxwell nanofluid. The Brownian motion produces random movement between fluid particles, which generate more heat, and as a result, the fluid temperature increases.

spinning disk with suction/injection effects has been examined. The flow is studied in the context of magnetization and radiation. The numerical results are found through bvp4c, while for comparison purposes, the computation is carried out via the RK4 technique. Figure 8a,b accordingly represents the action of the concentration field as a function of the thermophoresis term and Brownian. It can be shown that with *Nt*, the nanofluid concentration field raises while with *Nb*, the concentration profile drops. The suction

From the above studies the following conclusions have been drawn:

parameter (*Ws*) < 0 and Schmidt number *Sc* effect are illustrated in Figures 8c and 8d respectively. The concentration field reduces with an increase of both suction velocity and Schmidt number. The consequences of different parameters (*Bt*, *Ws*) on the axial and radial magnetic strength profiles are illustrated through Figure 9a–d respectively. It can observe that both parameters Batchlor number *Bt* and injection parameter *Ws* positively effect the magnetic strength profile along axial and radial direction.

Figure 10a,b express the nature of heat −*θ* 0 (0) and mass transfer −*φ* 0 (0) against radiation parameter *Rd* and Brownian motion parameter *Nb*, respectively. Because of the improving effect of radiation, the fluid temperature also rises, which enhances the heat transmission rate. Table 1 displays the numerical outcomes for skin friction and compared it with the published literature. Table 2 illustrates the comparison of bvp4c and RK4 techniques for the numerical outcomes. The numerical outputs for Sherwood number *Sh<sup>r</sup>* and Nusselt number *Nu<sup>r</sup>* are plotted in Table 3.


**Table 1.** The numerical outcomes for skin friction *F* 0 (0).


**Table 3.** The comparison of *RK*4 and Bvp4c for Sherwood *Sh<sup>r</sup>* and Nusselt number *Nur*, while keeping *ω* = 1.0, *δ* = 0.5, *β* = 0.1, *M* = 1.1, *Nb* = 0.5, *Nt* = 0.7, *Q<sup>w</sup>* = 1.2 and *Rd* = 0.5.


#### **6. Conclusions**

In the present mathematical model, the Maxwell nanoliquid flow over a porous spinning disk with suction/injection effects has been examined. The flow is studied in the context of magnetization and radiation. The numerical results are found through bvp4c, while for comparison purposes, the computation is carried out via the RK4 technique. From the above studies the following conclusions have been drawn:


**Author Contributions:** M.B. wrote the original manuscript and performed the numerical simulations. M.A.K. reviewed the mathematical results and S.-S.Z., T.M. and M.A.K. restructured the manuscript and revised the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deanship of Scientific Research at King Khalid University, Abha, Saudi Arabia for funding this work through research groups program under grant number R.G.P-1/142/42 and the APC was funded by S.S.Z.

**Acknowledgments:** The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University, Abha, Saudi Arabia for funding this work through research groups program under grant number R.G.P-1/142/42.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

