**Numerical Scrutinization of Darcy-Forchheimer Relation in Convective Magnetohydrodynamic Nanofluid Flow Bounded by Nonlinear Stretching Surface in the Perspective of Heat and**

#### **Ghulam Rasool 1, \* , Anum Shafiq 2 , Marei S. Alqarni 3,4 , Abderrahim Wakif 5 , Ilyas Khan 6 and Muhammad Shoaib Bhutta 1**


#### **1. Introduction**

A simple base fluid, for example, water, ethylene glycol, and oil, etc., when upgraded with the suspension of nanometric metallic strong conductive particles, is termed as nanofluid. Such a formulation sufficiently intensifies the conduction abilities of the base fluid. Numerous applications have been discovered for the so-called nanofluids in the industrial and engineering aspects as well as in bio-medicine. For example, vehicle cooling, heat exchangers, cooling and transformer cooling, electronic cooling, and many others are typical and widely used applications of nanofluids. These are also applicable in the medical treatments, especially cancer and tumor treatments, resonance imaging, and wound treatment, etc., and they are typically dependent on the conductive nature of nanofluids. Choi [1] introduced the definition of nanofluid in his experimental work where he proved that the suspension of nanoparticles in typical fluids drastically changes the thermo-physical properties of the fluid. Later on, Buongiorno [2] modeled the same concept in the perspective of convective transport of nanofluids. Adding details to the concept of nanofluids, Buongiorno emphasized the fact that the Brownian diffusion and Thermophoresis are two major slip factors in the transport of nanofluids. Afterwards, several interesting attempts have been reported by renowned researchers of fluid mechanics. For instance, Khan et al. [3] reported convection phenomena in nanofluid flow passing a linear stretching surface using the Keller–Box numerical method for the final solutions of the modeled governing problems. An interesting study of Mustafa et al. [4] disclosed an analysis on stagnation spot flow of nanofluids involving linear stretched sheet. For more details on this topic, one can see [5–24] and cross references cited therein.

Flow past a linear as well as nonlinear stretching surface relates the fluid mechanics with several important industrial and engineering setups, such as hot rolling, polymeric extrusion, continuously stretching done in plastic thin films, crystal growth, fiber production, and metallic extrusion, etc. Numerous articles are available in the literature that explain the flow caused by stretching surfaces, whether linear or nonlinear rates. The importance of linear stretching cannot be neglected but the flow caused by non-linear stretching rates

**Citation:** Rasool, G.; Shafiq, A.; Alqarni, M.S.; Wakif, A.; Khan, I.; Bhutta, M.S. Numerical Scrutinization of Darcy-Forchheimer Relation in Convective Magnetohydrodynamic Nanofluid Flow Bounded by Nonlinear Stretching Surface in the Perspective of Heat and Mass Transfer. *Micromachines* **2021**, *12*, 374. https://doi.org/10.3390/mi12040374

Academic Editor: Junfeng Zhang and Ruijin Wang

Received: 9 March 2021 Accepted: 24 March 2021 Published: 1 April 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/).

have always played a significant role in the above mentioned procedures, especially in polymeric extrusion. In all such scenarios, the work of Cortell [25] is considered to be a pioneer. He considered a viscous fluid for studying heat and mass transfer developments driven by a nonlinear or linear stretching in the sheet. The prescribed wall temperature and constant wall temperature were both discussed in this study. Vajravelu [26] reported a study on exploration of heat transfer developments in a viscous fluid flow past a stretching sheet surface using power law velocity distribution with nonlinear stretching rate. Rana and Bhargava [27] reported a study on the fluid flow analysis and heat transfer aspects involving a nonlinear stretching rate in nanofluids flow over a sheet.

The type of nanofluids governed in the categorical classification of non-Newtonian fluids have gained special attention for numerous engineering and industrial applications based on their extended contributions in nuclear, chemical, metallic, polymeric, and plastic industries. Shampoos, paints, apple sauce, katchup, as well as different type of oils are typical genuine examples of non-Newtonian fluids. Having extended viscosity, the non-Newtonian fluids are best treated by involving the viscoelastic terms (second grade model), which is known as sub-categorical classification of differential non-Newtonian liquids related with the normal stress attribute. For a better understanding of this model, one can read [28–34] and the references cited therein.

The studies mentioned above are either concerned with a linear stretching surface with effective convective heating or nonlinear stretching sheets without the involvement of convective conditions together with the Darcy–Forchheimer model. Here, for the first time, we involve convective conditions to consider a visco-elastic and strictly incompressible nanoliquid (nanofluid) flow bounded by a nonlinear flat stretching surface. Firstly, the model shapes in mathematical form using the famous Navier stokes equations for incompressible non-Newtonian nanofluid. The leading problems are then transformed into highly nonlinear ordinary problems via suitable transformations. A numerical scheme is implemented for finding the final solutions. From now onward, fluid means means incompressible viscoelastic nanofluid. The next section will physically justify the existence of the problem and mathematical expressions with properly defined boundary conditions. The rest of the article comprises of the results and discussion, a graphical display, and concluding remarks.

#### **2. Problem Formulation**

In this numerical investigation, we have invoked the convective boundary on the flow of nanofluid passing over a nonlinear flate stretching surface. These conditions are invoked to balance the temperature difference within the system. The system relies on Darcy–Forchheimer medium saturated via nanofluid over the stretching surface. Cartesian coordinates are considered to analyze the fluid flow. The flow direction is assumed along positive *x*−direction, whereas no-movement is allowed towards vertical. An induced magnetic effect is directly invoked to the surface normal direction via MHD; however, a tiny Reynolds number helps to dismiss the magnetic impact. The nonlinear stretching rate is taken into account via *n* as a positive integer, whereas the stretching velocity is assumed as *u* = *U<sup>w</sup>* = *mx<sup>n</sup>* at the bottom line. However, the velocity diminishes away from the surface and attains zero value at free surface (*u* = 0). The coefficient of heat and mass flux (*h*<sup>1</sup> and *h*2) are involved in convective boundary conditions. Figure 1 sketches a physical display.

**Figure 1.** Physical model and coordinate system.

The modified Navier stokes problems are given below:

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

$$\begin{split} v\frac{\partial u}{\partial y} + u\frac{\partial u}{\partial x} &= \nu \frac{\partial^2 u}{\partial y^2} - \left[ \frac{B\_0^2 \sigma}{\rho\_f} + \frac{\nu}{K} \right] u - \frac{\mathsf{C}\_b}{K^{1/2}} u^2 \\ &- k\_1 \left( u \frac{\partial^3 u}{\partial x \partial y^2} + \frac{\partial u}{\partial x} \frac{\partial}{\partial y} \left( \frac{\partial u}{\partial y} \right) + v \frac{\partial}{\partial y} \left( \frac{\partial^2 u}{\partial y^2} \right) - \frac{\partial u}{\partial y} \frac{\partial}{\partial x} \left( \frac{\partial u}{\partial y} \right) \right), \end{split} \tag{2}$$

$$a\frac{\partial T}{\partial y} + u\frac{\partial T}{\partial x} = a\frac{\partial}{\partial y}\left(\frac{\partial T}{\partial y}\right) + \frac{(\rho c)\_p}{(\rho c)\_f} \left[D\_{\text{Br}}\left(\frac{\partial \mathcal{C}}{\partial y}\frac{\partial T}{\partial y}\right) + \frac{D\_{Th}}{T\_{\infty}}\left(\frac{\partial T}{\partial y}\right)^2\right] - \frac{1}{\rho c\_f}\frac{\partial q\_r}{\partial y},\tag{3}$$

where *qr* is known as radiative heat flux. Given by Rosseland's approximation:

$$q\_r = -\frac{\partial (T^4)}{\partial y} \left(\frac{4\sigma^\*}{3k^\*}\right)\_{\prime} \tag{4}$$

where *σ* ∗ =Stefan–Boltzmann constant and *k* ∗ = mean absorption, respectively. Here, in,

4*T* 3 <sup>∞</sup>*T* − 3*T* 4 ∞ ∼ = *T* 4 , (5)

confirms,

$$\frac{\partial q\_r}{\partial y} = -\frac{\partial^2 T}{\partial y^2} \left( \frac{16\sigma^\* T\_\infty^3}{3k^\*} \right). \tag{6}$$

Therefore, the energy equation has the following final form:

$$\ln \frac{\partial T}{\partial \mathbf{x}} + \upsilon \frac{\partial T}{\partial y} = \varkappa \frac{\partial}{\partial y} \left( \frac{\partial T}{\partial y} \right) + \frac{(\rho c)\_p}{(\rho c)\_f} \left[ \frac{D\_{Th}}{T\_{\infty}} \left( \frac{\partial T}{\partial y} \right)^2 + D\_{Br} \left( \frac{\partial \mathbb{C}}{\partial y} \frac{\partial T}{\partial y} \right) \right] + \frac{\partial^2 T}{\partial y^2} \left( \frac{16 \sigma^\* T\_{\infty}^3}{3k^\*} \right) \frac{1}{(\rho c)\_f}, \tag{7}$$

$$v\frac{\partial \mathbf{C}}{\partial y} + u\frac{\partial \mathbf{C}}{\partial x} = D\_{Br} \left(\frac{\partial}{\partial y}\frac{\partial \mathbf{C}}{\partial y}\right) + \frac{D\_{Th}}{T\_{\infty}} \left(\frac{\partial}{\partial y}\frac{\partial T}{\partial y}\right) \tag{8}$$

Equations (1),(2),(7), and (8) are known as the governing equations with the following boundary conditions, as per the present model,

$$\text{If } u = \text{U}\_{\overline{w}} = m\mathbf{x}^n \text{, } h\_1(T\_f - T) = -k \frac{\partial T}{\partial y'} \text{, } v = 0, \text{ } h\_2(\mathbb{C}\_f - \mathbb{C}) = -D\_{\text{Br}} \frac{\partial \mathbb{C}}{\partial y} \quad \text{at} \quad y = 0, \tag{9}$$

$$
\Delta u = 0, \,\frac{\partial u}{\partial y} = 0, \,\text{C} = \text{C}\_{\infty}, \,\, T = T\_{\infty} \quad \text{as} \quad y \to \infty. \tag{10}
$$

In the above governing equations, *µ* is taken as dynamic viscosity, *ν* = *µ*/*ρ<sup>f</sup>* is known as kinematic viscosity, *ρ<sup>f</sup>* is taken as density, *σ* is involved for electrical conductivity, *α* = *k*/(*ρc*)*<sup>f</sup>* is well known thermal diffusive force, *k* is typical thermal conductiveness, *h*<sup>1</sup> = *hpx n*−1 <sup>2</sup> and *h*<sup>2</sup> = *hqx n*−1 <sup>2</sup> are known as non-uniform heat and mass (transfer) coefficients, respectively, *τ* is known as the ratio given between the effective nanoparticles heating capacity versus effective fluid heating capacity, *DBr* is known for Brownian diffusion, *DTh* is thermophoretic factor, and *B*<sup>0</sup> is magnetic effect. Defining,

$$v = -\frac{1}{2} \mathbf{x}^{\frac{n-1}{2}} \left[ \left( \frac{n-1}{n+1} \right) \frac{\partial f}{\partial \eta} \eta + f(\eta) \right] \sqrt{2m\nu(n+1)}, \quad u = mx^n \frac{\partial f}{\partial \eta},$$

$$\phi(\eta) = \frac{\mathbb{C} - \mathbb{C}\_{\infty}}{\mathbb{C}\_{w} - \mathbb{C}\_{\infty}}, \quad \theta(\eta) = \frac{T - T\_{\infty}}{T\_{w} - T\_{\infty}}, \quad \eta = \frac{1}{2} \sqrt{\frac{2\rho\_{ff}m(n+1)}{\mu}} x^{\frac{n-1}{2}} y. \tag{11}$$

Using (15) in (1), (2), (9) and (10) we have

$$\begin{aligned} &f'''' + ff'' + k\_0 \left[ \left(\frac{n+1}{2n}\right) f' f^{iv} + \left(f''\right)^2 - \left(\frac{n+1}{2n}\right) f'^{iv} f' \right] \\ &- \left(\frac{2n}{1+n}\right) (1+F\_r) \left(f'\right)^2 - \lambda \left(\frac{2}{n+1}\right) f' - M^2 \left(\frac{2}{n+1}\right) f' = 0, \end{aligned} \tag{12}$$

$$(1 + 4/3Rd)\theta'' + Prf\theta' + NbPr\theta'\phi' + NtPr\left(\theta'\right)^2 = 0,\tag{13}$$

$$
\phi'' + LePrf\phi' + \frac{Nt}{Nb}\theta'' = 0,
\tag{14}
$$

$$\begin{aligned} f(0) &= 0, & f' &= 1, & \phi'(0) &= -\gamma\_2(1 - \phi(0)), & \theta'(0) &= -\gamma\_1(1 - \theta(0)), \\ f'(\infty) &= 0, & f''(\infty) &= 0, & \phi(\infty) &= 0, & \theta(\infty) &= 0. \end{aligned} \tag{15}$$

Here, *γ<sup>i</sup>* for *i* = 1, 2 are the convective parameters extracted from *h*<sup>1</sup> = *hpx n*−1 <sup>2</sup> and *h*<sup>2</sup> = *hqx n*−1 <sup>2</sup> known as non-uniform heat transfer coefficient and mass transfer coefficient, respectively. *F<sup>r</sup>* is used for inertial frame, *M* is involved for magnetic impact, *λ* is treated as porosity, *Pr* is the typical Prandtl, *Le* is the well known Lewis number, *k*<sup>0</sup> is taken as viscoelastic parameter, and *Nb* and *Nt* Brownian diffusion and thermophoretic factors, respectively. Mathematically,

$$\gamma\_1 = \frac{h\_p}{k} \sqrt{\frac{2\nu}{m(n+1)}}, \gamma\_2 = \frac{h\_q}{D\_{Br}} \sqrt{\frac{2\nu}{m(n+1)}}, \quad F\_r = \frac{\mathbf{C}\_b \mathbf{x}}{K^{1/2}}, \quad M^2 = \frac{2\sigma B\_0^2}{m\mathbf{x}^{n-1}\rho\_f},$$

$$\lambda = \frac{2\nu}{Km\mathbf{x}^{n-1}}, \quad \text{Pr} = \frac{\nu}{\mathbf{a}'}, \quad L\mathbf{e} = \frac{a}{D\_{Br}}, \quad k\_1 = \frac{k\_0 m \mathbf{x}^{n-1}}{\nu}, \tag{16}$$

$$N\mathbf{b} = \frac{\left(\rho\_f c\right)\_p D\_{Br} (\mathbf{C}\_w - \mathbf{C}\_{\infty})}{\left(\rho\_f c\right)\_f \nu}, \quad N\mathbf{t} = \frac{\left(\rho\_f c\right)\_p D\_{Th} (T\_w - T\_{\infty})}{\left(\rho\_f c\right)\_f \nu T\_{\infty}}.$$

Furthermore, using *Re<sup>x</sup>* = *mxn*+1/*ν*, the local tiny Reynolds, the non-dimensional forms of physical quantities, are given below:

$$\begin{aligned} (\text{Re}\_x)^{1/2} \left( \frac{n+1}{2} \right)^{-1/2} \mathbb{C}\_f &= [(1-3\mathbf{k}\_0)f''(0)], \\ (\text{Re}\_x)^{-1/2} \left( \frac{n+1}{2} \right)^{-1/2} \mathrm{Nu} &= -\left(1+\frac{4}{3}\right) \mathrm{Rd} \left[ \theta'(0) \right], \\ (\text{Re}\_x)^{-1/2} \left( \frac{n+1}{2} \right)^{-1/2} \mathrm{Sh} &= -\left[\phi'(0) \right], \end{aligned} \tag{17}$$

#### **3. Methodology**

RK45 is one of the most frequently used numerical methods for solving IVPs for its accuracy and efficiency. The built-in codes involve basic concept of converting Boundary value problems into the initial value problems and subsequent problems are solved using a parallel scheme with RK45, such as shooting technique, secant method, etc., using an appropriate set of initial guesses to approximate the solutions. Herein, the boundary value problems are converted into initial value problems together with the given boundary conditions and thereafter, careful selection of initial guess for repeated iterations of the numerical scheme are chosen to obtain the solutions.

#### **4. Results and Discussion**

In this section, we have described the consequences noticed in flow profiles for all relevant parameters that are involved in leading problems of present nanofluid model. The arguments are comprehended by physical justifications as to why the variation is occurred and what consequence is being witnessed. The results are obtained using the numerical method and the data are graphically plotted in increasing and decreasing curves to witness the difference. The first section of the graphs belongs to the velocity profile, the second to temperature distribution, and third is related with the concentration of nanoparticles. In particular, Figure 2 is plotted to see the variation in velocity profile against augmented trend of visco-elastic parameter *k*0. The physical reasoning of the behavior noted for this parameter is connected with the combination of *k*<sup>1</sup> and kinematic viscosity *µ* having an inverse relation with fluid viscosity. For larger values of *k*0, the fluid viscosity shows a decline, and a consequent increasing trend is witnessed in the relative profile. The magnetic parameter *M* that is given in Figure 3 with a relevant variation in fluid motion produces a reduction. Physically, the affect of magnetic field is always in anti-directional with the fluid flow because of the surface normal bumps. The certain effect of magnetic field in the direction normal to the fluid flow is, for sure, the reason behind this trend of decline in fluid velocity. The similar trend in velocity profile can be figured out for the Forchheimer number *F<sup>r</sup>* that is given in Figure 4, but this time the decline is related with a higher friction and retardation offered by porous medium. The larger the Forchheimer number, the higher the frictional force that is offered to fluid motion and the consequence is a decline in fluid velocity. The thermal layer receives significant modification in the incremental direction for augmented values of thermal radiation *R<sup>d</sup>* given in Figure 5. The influence of thermal radiation eases the way of heat flux with a more convenient convection and the corresponding profile receives significant enhancement. Physically, thermal radiation is responsible for the enhancement of the thermal state of the fluid. Elevation in thermal Biot number *γ*<sup>1</sup> apparently results in a significant increase of the thermal layer and associated boundary layer thickness enhances, as shown in Figure 6. The heat flux coefficient that is involved in constitution of the said parameter is a sufficient justification for this behavior of thermal profile. Thermophoresis and Brownian motion are the two important factors of nanofluid flow and boundary layer phenomenon. By definition, Brownian motion is an uncertain movement of some particles (nanoparticles in this case) in the given medium for any time *t* > 0. The parameters are both

closed influential to each other as well as on the flow profiles. Here, the thermophoresis enforces a stronger thermophoretic push to the given nanoparticles and the nanoparticles receive certain inpredictive enhanced motion and, therefore, a disturbance appears and thermal profiles increase for both of the given factors shown in Figures 7 and 8, respectively. Figures 9–12 are related with the third important flow profiles, the concentration profiles, which is solely based on the concentration of nanoparticles. In particular, Figure 9 presents the display of variation that is noticed in concentration profile with respect to the elevation in Schmidt number. Physically, the inverse relationship between Brownian diffusion and kinematic viscosity is sufficient to justify the reducing trend noticed in concentration profile for augmented values of Schmidt number. The solute Biot number is an enhancing factor the concentration profile that is given in Figure 10. The coefficient of mass flux involved in the constitutive expression justifies this variation. An enhancement is noticed for elevated values of thermophoresis given in Figure 11. A decline can be seen in concentration profile for augmented Brownian diffusion parameter given in Figure 12. For sure, the decline is based on the in-predictive and uncertain fluctuation of nanoparticles. Figures 13 and 14 are the contour graphs sketched for both the linear and non-linear case, respectively, while Figures 15 and 16 are the density graphs. A minor but significant difference can be seen in both the linear and non-linear case settled at *n* = 1 & 1.5, respectively. Relevant data for skin-friction, heat, and mass flux numbers are given in Tables 1 and 2, respectively. Skin friction jumps for larger porosity and a larger Forchheimer number, therefore, declines the easy movement of fluid.The heat flux and mass flux both receive a reduction for the augmented values of Forchheimer number. Thermal radiation is an enhancing factor for both of the flux rates. An opposite trend is noticed in both Nusselt and Sherwood numbers for thermal Biot number. Heat flux enhances, while mass flux reduces with the passage of time.

**Figure 2.** Visco-elastic parameter versus velocity field.

**Figure 4.** Forchheimer number versus velocity field.

**Figure 5.** Thermal radiation parameter versus thermal profile.

**Figure 6.** Thermal Biot number versus thermal profile.

**Figure 7.** Thermophoresis versus thermal profile.

**Figure 8.** Brownian diffusion versus thermal profile.

**Figure 9.** Schmidt number versus concentration profile.

**Figure 10.** Solute Biot number versus concentration distribution.

**Figure 11.** Thermophoresis versus concentration distribution.

**Figure 12.** Brownian diffusion versus concentration distribution.

For Linear Case

**Figure 13.** Contour graph for linear case.

For non-Linear Case

**Figure 14.** Contour graph for nonlinear case.

#### Density Plot for Linear Case Density Plot for non-Linear Case



**Author Contributions:** Conceptualization, G.R.; methodology, G.R., A.S.; software, G.R., A.S., A.W.; validation, G.R., A.S., A.W., M.S.B.; formal analysis, G.R., A.S., A.W., I.K., M.S.A., M.S.B.; investigation, G.R., A.S., A.W., I.K., M.S.B.; resources, M.S.A.; data curation, G.R., A.S., A.W., M.S.B.; writing—original draft preparation, G.R., A.S., A.W., I.K., M.S.B., M.S.A.; writing—review and editing, G.R., A.S., A.W., I.K., M.S.B., M.S.A.; visualization, G.R., A.S., A.W.,; supervision, M.S.A., I.K.; project administration, M.S.A.; funding acquisition, M.S.A. All authors have read and agreed to the published version of the manuscript.

**Data Availability Statement:** All data are available within this manuscript.

**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/128/42.

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

#### **Nomenclature:**


#### **References**


*Micromachines* **2021**, *12*, 374

*Micromachines* **2021**, *12*, 374
