*Article* **Significant Involvement of Double Diffusion Theories on Viscoelastic Fluid Comprising Variable Thermophysical Properties**

**Muhammad Sohail 1,\* , Umar Nazir <sup>1</sup> , Omar Bazighifan <sup>2</sup> , Rami Ahmad El-Nabulsi 3,4,5,\*, Mahmoud M. Selim 6,7 , Hussam Alrabaiah 8,9 and Phatiphat Thounthong <sup>10</sup>**


**Abstract:** This report examines the heat and mass transfer in three-dimensional second grade non-Newtonian fluid in the presence of a variable magnetic field. Heat transfer is presented with the involvement of thermal relaxation time and variable thermal conductivity. The generalized theory for mass flux with variable mass diffusion coefficient is considered in the transport of species. The conservation laws are modeled in simplified form via boundary layer theory which results as a system of coupled non-linear partial differential equations. Group similarity analysis is engaged for the conversion of derived conservation laws in the form of highly non-linear ordinary differential equations. The solution is obtained vial optimal homotopy procedure (OHP). The convergence of the scheme is shown through error analysis. The obtained solution is displayed through graphs and tables for different influential parameters.

**Keywords:** viscoelastic material; group similarity analysis; thermal relaxation time; parametric investigation; variable magnetic field; error analysis

## **1. Introduction**

Fluid flows over stretched surfaces have applications in several fields and has significant involvement in the practical usage of several items. Scientists and engineers have made efforts to explore their features and usage in different processes. The mathematical relations of non-Newtonian materials are different as compared with Newtonian material. These materials are divided into different categories according to their properties. An important non-Newtonian fluid is a second grade fluid [1–7]. It has the following constitutive relation:

$$
\pi^\* = -PI + \mu F\_1 + \beta\_1 F\_2 + \beta\_2 F\_1 \ast F\_1, \; \beta\_1 \ge 0, \; \mu \ge 0, \; \beta\_1 + \beta\_2 = 0.
$$

**Citation:** Sohail, M.; Nazir, U.; Bazighifan, O.; El-Nabulsi, R.A.; Selim, M.M.; Alrabaiah, H.; Thounthong, P. Significant Involvement of Double Diffusion Theories on Viscoelastic Fluid Comprising Variable Thermophysical Properties. *Micromachines* **2021**, *12*, 951. https://doi.org/10.3390/ mi12080951

Academic Editors: Lanju Mei and Shizhi Qian

Received: 21 June 2021 Accepted: 27 July 2021 Published: 12 August 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/).

Hayat et al. [1] studied the mixed convection in the second grade model over a stretching cylinder. They modeled the problem in two dimensions with thermal transport by taking the variable thermal conductivity. They used a homotopy method for the solution. They studied the contribution of several emerging parameters on the solution through graphs. They noticed the decrease in velocity field for a mixed convection parameter. Massoudi et al. [2] reported the study on the second grade model with temperature dependent viscosity between parallel plates. They presented a comparative study for the validity of obtained solution through tabular data. An exact solution through an oscillating sphere for a second grade model was computed by Fetecau et al. [3]. They found that the solution is periodic and it is independent of initial data. Hankel and Laplace transforms were engaged by Kamran et al. [4] to handle the modeled equations for fractional second grade model in cylindrical coordinates. They presented that the fractional model present the fluid flow phenomenon more accurately as compared with the ordinary derivatives. Chauhan and Kumar [5] studied the unsteady mechanics for second grade model in partially filled porous channel. They used a Laplace transform technique to analyze the solution. They observed the increase in velocity field against time parameter. A rotating viscoelastic model with ramped wall temperature condition for exact solution was reported by Mohamad et al. [6]. They plotted the solution against numerous emerging parameters. They noticed the dual behavior of velocity against the rotation parameter. Hayat et al. [7] examined the comportment of chemical reaction with solutal and thermal transport in a second grade model passed over a bi-directional stretched surface. They found the increase in dimensionless stress against ratio parameter and viscoelastic parameter. Moon et al. [8] discussed the phenomenon of heat transfer and Weber number including droplets of xanthan gum solution (non-Newtonian) and DI-water (Newtonian) over a heated surface. They noted that the DI-water droplet has higher spreading diameter as compared with the non-Newtonian (droplet) because of variation in fluid difference. German and Bertola [9] studied free-fall related to the liquid drops due gravity of force. They imagined high speed of drops based on viscoelastic fluids. They found that shape of the drop is changed under the action of yield stress. An and Lee [10] experimentally discussed the oscillations (free falling) of drops (shear-thinning) due to the force of gravity based on viscoelastic fluids passed over a solid surface. Moon et al. [11] suggested a mixed regime (coalescence occurs) using sequential images and mixed regime is exaggerated due to volumes of static droplet volumes, static droplets and Weber numbers. They estimated the film thickness (between two drops) via lubrication theory. Zhao and Khayat [12] discussed the flow behavior in view of shear-thickening and shear-thinning of a jet (non-Newtonian) over a flat plate via a hydraulic jump. Moon et al. [13] considered the features of droplets (non-Newtonian) on solid surfaces considering various Weber numbers. They considered xanthan gum solution to produce the droplets (non-Newtonian) measured via spreading diameters and camera (high speed).

Transportation of heat has applications in different engineering aspects and it is now a hot topic for researchers working in the field of engineering and applied mathematics. Several researchers are working on transport phenomena actively. For instance, Naseem et al. [14] worked on third grade nanofluid passed over a Riga plate. They considered several physical effects while modelling the transport equations and the resulting equations were simplified via boundary layer theory with the solution approximated analytically. They noticed the rise in velocity profile for modified magnetic parameter and Reynolds number. Sahoo and Poncet [15] studied the Blasius flow of a fourth grade model with heat transfer in a porous permeable stretching surface. Heat transfer in MHD viscous stagnation point dusty fluid with a non-uniform heat source in a porous stretching sheet was studied by Ramesh et al. [16]. They presented a comparative analysis and solved the resulting equations with the help of a shooting method. They discussed the contribution of several parameters on velocity and temperature fields. They noticed the decline in thermal field against higher Prandtl number. Qiu et al. [17] studied the thermal transport in channel via finite volume technique. They also reported the analysis of entropy in their

findings. They discussed the impact of nanoparticles volume fraction on flow and entropy. Khan et al. [18] studied the heat transfer in a stagnation point Powell–Eyring model in a stretching cylinder with variable properties analytically. They observed the enhancement in temperature and velocity fields against the growing values or curvature parameter. Also they listed the numerical values for heat transfer coefficient against different parameters. Few recent studies covering non-linear transport problems with different effects have been reported in [19–25].

The objective of current inspection is to analyze the comportment of variable properties in heat and mass transportation in a second grade steady incompressible model past over a bi-directional elongating surface. This article is organized as follows: Section 1 contains the literature survey; modeling of considered problem is included in Section 2 with important physical quantities; Section 3 comprises methodology and results with key findings covered in Sections 4 and 5 respectively.

#### **2. Mathematical Drafting of Viscoelastic Fluid with Thermal and Mass Transport**

An analysis of the transport phenomenon on a second grade fluid [7,14] over a bidirectional elastic surface is presented in Figure 1. It is assumed that the sheet is stretched along x- and y- directions, respectively, and flow occupies the region normal to *x*- and *y*-axis. The sheet is kept at temperature "*Tw*" and concentration "*Cw*". Along *x*-axis, the velocity is "*U<sup>W</sup>* = *ax*" and "*V<sup>W</sup>* = *by*" is along the y-axis. The following important considerations have been adopted to derive the conservation laws

**Figure 1.** Geometry of second grade fluid model.


The flowing resulting equations [7] appears by using the above stated assumptions

$$u\_x + v\_y + w\_z = 0,\tag{1}$$

$$\begin{split} \mu u\_x + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} - \theta \frac{\partial^2 u}{\partial z^2} \\ + a\_0 \left[ u \frac{\partial^3 u}{\partial x \partial z^2} + w \frac{\partial^3 u}{\partial z^3} - \frac{\partial u}{\partial x} \frac{\partial^2 u}{\partial z^2} - \frac{\partial u}{\partial z} \frac{\partial^2 w}{\partial z^2} - 2 \frac{\partial u}{\partial z} \frac{\partial^2 u}{\partial x \partial z} \\ - 2 \frac{\partial w}{\partial z} \frac{\partial^2 u}{\partial z^2} \right] + \frac{\sigma}{\rho} B\_a^2(x, y) u = 0, \end{split} \tag{2}$$

$$\begin{split} u v\_{\lambda} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} - \theta \frac{\partial^2 v}{\partial z^2} \\ + a\_0 \left[ v \frac{\partial^3 v}{\partial x \partial z^2} + w \frac{\partial^3 v}{\partial z^3} - \frac{\partial v}{\partial x} \frac{\partial^2 v}{\partial z^2} - \frac{\partial v}{\partial z} \frac{\partial^2 w}{\partial z^2} - 2 \frac{\partial v}{\partial z} \frac{\partial^2 v}{\partial x \partial z} \\ - 2 \frac{\partial w}{\partial z} \frac{\partial^2 v}{\partial z^2} \right] + \frac{\sigma}{\rho} B\_a^2(x, y) v = 0, \end{split} \tag{3}$$

$$\begin{aligned} \mu T\_{\mathbf{x}} + vT\_{\mathbf{y}} + wT\_{\mathbf{z}} + \mathfrak{a}\_{d} \begin{bmatrix} \left(\mu u\_{\mathbf{x}} + vu\_{\mathbf{y}} + wu\_{\mathbf{z}}\right)T\_{\mathbf{x}} + \left(\mu v\_{\mathbf{x}} + v\frac{\partial v}{\partial \mathbf{y}} + w\frac{\partial v}{\partial \mathbf{z}}\right)T\_{\mathbf{y}} \\ + \left(\mu w\_{\mathbf{x}} + v\frac{\partial w}{\partial \mathbf{y}} + w\frac{\partial w}{\partial \mathbf{z}}\right)T\_{\mathbf{z}} + 2\mu vT\_{\mathbf{x}\mathbf{y}} + 2\nu wT\_{\mathbf{y}\mathbf{z}} \\ + 2\mu wT\_{\mathbf{x}\mathbf{z}} + u^{2}T\_{\mathbf{x}\mathbf{x}} + v^{2}T\_{\mathbf{y}\mathbf{y}} + w^{2}T\_{\mathbf{z}\mathbf{z}} \end{bmatrix} \end{aligned} \tag{4}$$
 
$$\begin{aligned} -\nabla[\mathbf{K}\_{A}(T)\nabla\mathbf{T}] = \mathbf{0} \end{aligned} \tag{5}$$

$$\begin{aligned} \left(\boldsymbol{u}\mathbf{C}\_{x} + \boldsymbol{v}\mathbf{C}\_{y} + \boldsymbol{w}\mathbf{C}\_{z} + \boldsymbol{a}\_{b} \begin{bmatrix} \left(\boldsymbol{u}\boldsymbol{u}\_{x} + \boldsymbol{v}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{y}} + \boldsymbol{w}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{z}}\right)\mathbf{C}\_{x} + \left(\boldsymbol{u}\boldsymbol{v}\_{x} + \boldsymbol{v}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{y}} + \boldsymbol{w}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{z}}\right)\mathbf{C}\_{y} \\ + \left(\boldsymbol{u}\boldsymbol{w}\_{x} + \boldsymbol{v}\frac{\partial\boldsymbol{w}}{\partial\boldsymbol{y}} + \boldsymbol{w}\frac{\partial\boldsymbol{w}}{\partial\boldsymbol{z}}\right)\mathbf{C}\_{z} + 2\boldsymbol{u}\boldsymbol{v}\mathbf{C}\_{xy} + 2\boldsymbol{v}\boldsymbol{w}\mathbf{C}\_{yz} \\ + 2\boldsymbol{u}\boldsymbol{w}\mathbf{C}\_{xz} + \boldsymbol{u}^{2}\mathbf{C}\_{xx} + \boldsymbol{v}^{2}\mathbf{C}\_{yy} + \boldsymbol{w}^{2}\mathbf{C}\_{zz} \end{bmatrix} \end{aligned} \tag{5}$$
 
$$\begin{aligned} -\nabla\left[\mathbf{D}\_{\mathbf{A}}(\mathbf{T})\nabla\mathbf{C}\right] = 0. \end{aligned} \tag{6}$$

Boundary conditions for the dimensional problem are

$$\begin{cases} \text{ } \boldsymbol{u} = \mathcal{U}\_{\mathcal{W}} = a\mathbf{x}, \; \boldsymbol{v} = V\_{\mathcal{W}} = by\_{\prime}\boldsymbol{w} = \mathbf{0}, \; T = T\_{\mathcal{W}}. \; \mathbb{C} = \mathbb{C}\_{\mathcal{W}} \text{ at } \boldsymbol{z} = \mathbf{0}. \\\ \boldsymbol{u} \to \mathbf{0}, \; \boldsymbol{v} \to \mathbf{0}, \; T \to T\_{\infty \prime} \; \mathbb{C} \to \mathbb{C}\_{\infty} \; for \; \boldsymbol{z} \to \infty. \end{cases} \tag{6}$$

With the use of the following similarity variables, the governing law reduce to

$$\begin{cases} \boldsymbol{u} = a \boldsymbol{x} f'[\eta], \; \boldsymbol{v} = a \boldsymbol{y} g'[\eta], \; \boldsymbol{w} = -(a\theta)^{\frac{1}{2}} [f[\eta] + \boldsymbol{g}[\eta]], \\\ \boldsymbol{\theta}[\eta] = \frac{\boldsymbol{T} - \boldsymbol{T}\_{\infty}}{\boldsymbol{T}\_{\mathrm{w}} - \boldsymbol{T}\_{\infty}}, \; \boldsymbol{\phi}[\eta] = \frac{\boldsymbol{\mathcal{C}} - \boldsymbol{\mathcal{C}}\_{\infty}}{\boldsymbol{\mathcal{C}}\_{\mathrm{w}} - \boldsymbol{\mathcal{C}}\_{\infty}}, \; [\eta] = \left(\frac{\boldsymbol{a}}{\theta}\right)^{\frac{1}{2}} \boldsymbol{z}, \end{cases} \tag{7}$$

$$\begin{aligned} -Mf'[\eta] - f'[\eta]^2 + (f[\eta] + g[\eta])f''[\eta] + f^{(3)}[\eta] + \mathcal{R}(f''[\eta](f''[\eta] - g''[\eta])) \\ -2(f'[\eta] + g'[\eta])f^{(3)}[\eta] + (f[\eta] + g[\eta])f^{(4)}[\eta]) = 0, \end{aligned} \tag{8}$$

$$\begin{aligned} -M\mathbf{g'}[\eta] - \mathbf{g'}[\eta]^2 + (f[\eta] + \mathbf{g}[\eta])\mathbf{g''}[\eta] + \mathbf{g}^{(3)}[\eta] + \mathcal{R}((f''[\eta] - \mathbf{g''}[\eta])\mathbf{g''}[\eta] \\ -2(f'[\eta] + \mathbf{g'}[\eta])\mathbf{g^{(3)}}[\eta] + (f[\eta] + \mathbf{g}[\eta])\mathbf{g^{(4)}}[\eta] &= 0, \end{aligned} \tag{9}$$

$$\begin{split} (f[\eta] + \mathfrak{g}[\eta])\theta'[\eta] + \frac{1}{\mathcal{P}r}(1 + \gamma\_1\theta[\eta])\theta''[\eta] - \mathfrak{a}\_1(f[\eta] + \mathfrak{g}[\eta])(f'[\eta] \\ + \mathcal{g}'[\eta]\theta'[\eta] + (f[\eta] + \mathfrak{g}[\eta])\theta''[\eta]) &= 0, \end{split} \tag{10}$$

$$\begin{aligned} (f[\eta] + \mathfrak{g}[\eta])\phi'[\eta] + \frac{1}{\mathfrak{S}c}(1 + \gamma\_2\mathfrak{e}[\eta])\phi''[\eta] - \mathfrak{a}\_2(f[\eta] + \mathfrak{g}[\eta])(f'[\eta] \\ + \mathfrak{g}'[\eta]\phi'[\eta] + (f[\eta] + \mathfrak{g}[\eta])\phi''[\eta]) &= 0, \end{aligned} \tag{11}$$

$$\begin{cases} \begin{array}{ll} f(0) = 0, \ g(0) = 0, \ f'(0) = 1, \ g'(0) = \delta, \ \theta(0) = 1, \ \phi(0) = 1, \\\ f'(\infty) = 0, \ g'(\infty) = 0, \ \theta(\infty) = 0, \ \phi(\infty) = 0. \end{cases} \end{cases} \tag{12}$$

#### *Physical Quantities*

The study of heat, mass transfer rates and dimensionless stress at the boundary has significant applications and usage in industry. Therefore, scientists and engineers are keenly observing their features against different physical parameters which influence them directly. These quantities are defined as:

$$\mathcal{C}\_{XF} = \frac{\left. \tau\_{\mathbf{x}\mathbf{z}} \right|\_{z} = 0}{\rho \left( u\_{\mathbf{w}} \right)^{2}}, \mathcal{C}\_{YF} = \frac{\left. \tau\_{\mathbf{y}\mathbf{z}} \right|\_{z} = 0}{\rho \left( v\_{\mathbf{w}} \right)^{2}}, \tag{13}$$

$$\tau\_{\mathbf{x}\mathbf{z}} = \left[\mu\frac{\partial\mu}{\partial z} + \mathfrak{a}\_0 \left[u\frac{\partial^2 u}{\partial x \partial z} + v\frac{\partial^2 u}{\partial y \partial z} + w\frac{\partial^2 u}{\partial z^2} + \frac{\partial u}{\partial x}\frac{\partial u}{\partial z} + \frac{\partial v}{\partial z}\frac{\partial v}{\partial x} - \frac{\partial w}{\partial z}\frac{\partial u}{\partial y}\right]\right]\_{z=0} \tag{14}$$

$$\pi\_{yz} = \left[ \mu \frac{\partial v}{\partial z} + \mu\_0 \left[ u \frac{\partial^2 v}{\partial x \partial z} + v \frac{\partial^2 v}{\partial y \partial z} + w \frac{\partial^2 v}{\partial z^2} + \frac{\partial u}{\partial y} \frac{\partial u}{\partial z} + \frac{\partial v}{\partial z} \frac{\partial v}{\partial y} - \frac{\partial w}{\partial z} \frac{\partial v}{\partial y} \right] \right]\_{z=0} \tag{15}$$

$$Nu\_{xy} = \frac{(\varkappa + y)Q\_w^\*}{K(T)[T\_w - T\_\infty]}, \ Q\_w^\* = -K(T)\frac{\partial T}{\partial z}|\_{z=0\prime} \tag{16}$$

$$\mathcal{S}u\_{xy} = \frac{(\mathbf{x} + \mathbf{y})M\_w^\*}{D\_B(T)[\mathbb{C}\_w - \mathbb{C}\_\infty]}, \ \mathcal{M}\_w^\* = -D\_B(T)\frac{\partial \mathcal{C}}{\partial z}|\_{z=0} \tag{17}$$

After boundary layer theory, the dimensionless form is:

$$\mathbf{C}\_{\rm Fx}^{\*} = f''(0) - \mathcal{R}[f(0) + \mathbf{g}(0)]f'''(0) + \mathcal{R}[f'(0) + \mathbf{g}'(0)]\left[f'(0) + 2\mathcal{R}f'(0)f''(0), \tag{18}$$

$$\mathbf{f}\_{\mathbf{F}y}^{\*} = \mathbf{g}''(0) - \mathbb{R}[f(0) + \mathbf{g}(0)]\mathbf{g}'''(0) + \mathbb{R}\left[f'(0) + \mathbf{g}'(0)\right]\mathbf{g}''(0) + 2\mathbb{R}f'(0)\mathbf{g}'(0)\tag{19}$$

$$\boldsymbol{H}\_{\rm xy}^{\*} = -\left(\boldsymbol{R}\_{\varepsilon \rm xy}\right)^{\frac{1}{2}} \boldsymbol{\theta}'(\mathbf{0}) \, \boldsymbol{M}\_{\rm xy}^{\*} = -\left(\boldsymbol{R}\_{\varepsilon \rm xy}\right)^{\frac{1}{2}} \boldsymbol{\phi}'(\mathbf{0}).\tag{20}$$

#### **3. Numerical Method for Solution**

*C* ∗

Modelling of the fluid flow problems results in the form of a set of coupled non-linear differential equations. The derived problem is highly non-linear and coupled. Due to high non-linearity, an exact solution is not possible. Researchers proposed several schemes to handle the non-linear complex differential equations. Here, optimal homotopy analysis procedure (OHAP) [7,14,19–23] is engaged due to its several advantages.

This section covers the necessary steps for the adopted procedure. It has the following steps:


The linear operators with initial guesses are:

$$\begin{cases} \begin{array}{c} L\_f^\* = \frac{D^3}{D\eta^3} - \frac{D}{D\eta}, \ L\_g^\* = \frac{D^3}{D\eta^3} - \frac{D}{D\eta}, \\\ L\_t^\* = \frac{D^2}{D\eta^2} - 1, \ L\_c^\* = \frac{D^2}{D\eta^2} - 1, \end{array} \tag{21}$$

$$\begin{cases} f\_q(\eta) = 1 - e^{-\eta} \text{ } \text{g}\_{\emptyset} = \gamma [1 - e^{-\eta}] \text{ } \\\ \theta\_q(\eta) = e^{-\eta} \text{ } \phi\_q(\eta) = e^{-\eta} \text{ } \end{cases} \tag{22}$$

The operators in Equation (21) obeys:

$$\begin{aligned} L\_f^\*[Q\_1 + Q\_2e^{-\eta} + Q\_3e^{\eta}] &= 0, \\ L\_g^\*[Q\_4 + Q\_5e^{-\eta} + Q\_6e^{\eta}] &= 0, \\ L\_t^\*[Q\_7e^{-\eta} + Q\_8e^{\eta}] &= 0, \\ L\_c^\*[Q\_9e^{\eta} + Q\_{10}e^{-\eta}] &= 0. \end{aligned} \tag{23}$$

where *Qb*(*b* = 1, 2, . . . , 10) are unknowns.

Using the concepts of minimization of average squared residual error [7–14,19–23]:

$$\delta\_m^f = \frac{1}{B+1} \sum\_{r=0}^B \left[ \mathbb{S}\_f \left( \sum\_{L=0}^a f(\eta)\_r \sum\_{L=0}^a \mathbb{S}(\eta) \right) \right]^2. \tag{24}$$

$$\delta\_m^{\mathcal{S}} = \frac{1}{B+1} \sum\_{r=0}^{B} \left[ \mathcal{S}\_{\mathcal{S}} \left( \sum\_{L=0}^{a} \mathcal{f}(\eta)\_{\prime} \sum\_{L=0}^{a} \mathcal{g}(\eta) \right) \right]^2. \tag{25}$$

$$\delta\_m^{\theta} = \frac{1}{B+1} \sum\_{r=0}^{B} \left[ \mathbb{S}\_{\theta} \left( \sum\_{l=0}^{a} \hat{f}(\eta)\_{l} \sum\_{L=0}^{a} \hat{g}(\eta)\_{l} \sum\_{L=0}^{a} \hat{\theta}(\eta) \right) \right]^2,\tag{26}$$

$$\delta\_m^\Phi = \frac{1}{B+1} \sum\_{r=0}^B \left[ S\_\Phi \left( \sum\_{L=0}^a f(\eta), \sum\_{L=0}^a \xi(\eta), \sum\_{L=0}^a \hat{\theta}(\eta), \sum\_{L=0}^a \hat{\phi}(\eta) \right) \right]^2,\tag{27}$$

where

$$
\delta\_i^t = \delta\_i^f + \delta\_i^g + \delta\_i^\theta + \delta\_i^\phi. \tag{28}
$$

The minimum error at second order is 0.00031589832079710834 and optimal values at third order are *B<sup>f</sup>* = −1.2160, *B<sup>g</sup>* = −1.1638, *B<sup>θ</sup>* = −0.9509, *B<sup>φ</sup>* = −0.5871, by fixing the involved parameters as *R* = 0.1, *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *M* = 0.1, *δ* = 0.8.

$$\begin{array}{l} f = 1.0 - 1.0e^{-z} - 0.6080M^2 + 0.6080e^{-z}M^2 - 0.6080R + 0.6080e^{-z}R \\ \qquad + 0.6080e^{-z}M^2z + 0.6080e^{-z}Rz - 0.4053\delta \\ \qquad + 0.2026e^{-2z}\delta + 0.2026e^{-z}\delta - 1.0133R\delta \\ \qquad - 0.4053e^{-2z}R\delta + 1.4187e^{-z}R\delta + 0.6080e^{-z}z\delta \\ \qquad + 0.6080e^{-z}Rz\delta \end{array} \tag{29}$$

$$\begin{array}{l} \rm g = 1.1939\delta + 0.1939e^{-2z}\delta - 1.3879e^{-z}\delta - 0.5819M^2\delta + 0.5819e^{-z}M^2\delta \\ \quad - 0.5819R\delta + 0.5819e^{-z}R\delta + 0.5819e^{-z}M^2z\delta \\ \quad + 0.5819e^{-z}Rz\delta - 0.5819\delta^2 + 0.5819e^{-z}\delta^2 - 0.9698R\delta^2 \\ \quad - 0.3879e^{-2z}R\delta^2 + 1.3578e^{-z}R\delta^2 + 0.5819e^{-z}z\delta^2 \\ \quad + 0.5819e^{-z}Rz\delta^2, \end{array} \tag{30}$$

*θ* = −0.3169*e* <sup>−</sup>2*<sup>z</sup>* + 1.3169*e* <sup>−</sup>*<sup>z</sup>* <sup>−</sup> 0.4754*<sup>e</sup>* −*z z* + 0.47543*<sup>e</sup>* −*z z Pr* − 0.3169*e* −2*z δ* + 0.3169*e* −*z δ* − 0.4754*e* −*z zδ* + 0.11887*e* <sup>−</sup>3*zα*<sup>1</sup> − 0.95097*e* <sup>−</sup>2*zα*<sup>1</sup> + 0.8321*e* <sup>−</sup>*zα*<sup>1</sup> <sup>−</sup> 0.9509*<sup>e</sup>* −*z zα*<sup>1</sup> + 0.3566*e* −3*z δα*<sup>1</sup> − 1.9019*e* −2*z δα*<sup>1</sup> + 1.5453*e* −*z δα*<sup>1</sup> − 1.4264*e* −*z zδα*<sup>1</sup> + 0.2377*e* −3*z δ* <sup>2</sup>*α*<sup>1</sup> <sup>−</sup> 0.9509*<sup>e</sup>* −2*z δ* <sup>2</sup>*α*<sup>1</sup> + 0.7132*e* −*z δ* <sup>2</sup>*α*<sup>1</sup> <sup>−</sup> 0.475497*<sup>e</sup>* −*z zδ* <sup>2</sup>*α*<sup>1</sup> <sup>−</sup> 0.3169*e* <sup>−</sup>2*zγ*<sup>1</sup> *Pr* + 0.3169*e* <sup>−</sup>*zγ*<sup>1</sup> , *Pr* (31) = −0.1957*e* <sup>−</sup>2*<sup>z</sup>* + 1.1957*e* <sup>−</sup>*<sup>z</sup>* <sup>−</sup> 0.2935*<sup>e</sup>* −*z z* + 0.2935*<sup>e</sup>* −*z z* Sc − 0.1957*e* −2*z δ* + 0.1957*e* −*z δ* − 0.2935*e* −*z zδ* + 0.0733*e* <sup>−</sup>3*zα*<sup>2</sup> − 0.5871*e* <sup>−</sup>2*zα*<sup>2</sup> + 0.5137*e* <sup>−</sup>*zα*<sup>2</sup> <sup>−</sup> 0.5871*<sup>e</sup>* −*z zα*<sup>2</sup> + 0.2201*e* −3*z δα*<sup>2</sup> − 1.1742*e* −2*z δα*<sup>2</sup> + 0.9540*e* −*z δα*2 − 0.8806*e* −*z zδα*<sup>2</sup> + 0.1467*e* −3*z δ* <sup>2</sup>*α*<sup>2</sup> <sup>−</sup> 0.5871*<sup>e</sup>* −2*z δ* <sup>2</sup>*α*<sup>2</sup> + 0.4403*e* −*z δ* <sup>2</sup>*α*<sup>2</sup> <sup>−</sup> 0.2935*<sup>e</sup>* −*z zδ* <sup>2</sup>*α*<sup>2</sup> <sup>−</sup> 0.1957*e* <sup>−</sup>2*zγ*<sup>2</sup> Sc + 0.1957*e* <sup>−</sup>*zγ*<sup>2</sup> Sc . (32)

## **4. Analysis and Discussion**

The assessments of vital applications of thermal energy and mass transport (using second grade liquid) in industrial and engineering areas are addressed, including various features. In this current problem, non-Fourier's theory is investigated in energy and mass transport equations along with the concept of variable properties (mass diffusion and thermal conductivity). The motion in nanoparticles is induced because of movement of a melting 3D-surface. The simulations of temperature, diffusion of mass and flow behavior are captured in the form of graphs and tables via an analytical approach. The error analysis is presented with the help of Table 1. The graphical discussions are addressed below:


**Table 1.** Computation of averaged squared residuals errors of velocity, temperature, and concentration solution.

**Assessments of flow phenomena via physical parameters:** the distribution of flow phenomena is analyzed with respect to magnetic number (*M*), second grade fluid number (*R*) and velocity stretching ratio number (*δ*). In this current problem, *M* is considered as variable magnetic number for measurement of flow behavior in both directions (vertical and horizontal). The variation strength of a magnetic field is considered during the flow of particles at the surface of the melting sheet. This effect is analyzed by Figures 2 and 3. The decreasing character in the motion of particles is noticed via enhancement strength regarding magnetic field. The declination in flow is due to a negative force called Lorentz force appeared in momentum equations. The retardation forces play a reducing role in motion of fluid particles. Therefore, a magnetic number is used to adjust (MBL) momentum boundary layer thickness. It is noticed that the last terms of dimensionless momentum equations represent negative Lorentz forces. These negative Lorentz forces generate hindrance during the flow of fluid particles. The reduction is investigated in velocity profiles for *M* = 0.0, 0.2, 0.4, 0.6, 0.8. Therefore, a local minimum trend is noticed for *M* = 0.0, 0.2, 0.4, 0.6, 0.8. The parameter related to *R* is known as a second grade number whereas the character of *R* is simulated on the flow behavior. In this case, the reduction in flow mechanism is conducted through the numerical values of *R* and this graphical impact is considered in Figures 4 and 5. It is noticed that the appearance of *R* is modeled due to the appearance of second grade fluid in the current model. Meanwhile, flow in the horizontal and vertical directions is reduced, taking higher values of *R*. Moreover, flow in the absence of *R* is higher as compared with flow situation in the presence of *R*. Simply, it is included that flow for the case Newtonian liquid becomes higher than as compared flow for the case of non-Newtonian liquid. Figure 6 provides the role of *δ* on distribution of velocities. The enhancement in velocity profiles is conducted using higher values of *δ* while the parameter related to *δ* is addressed as velocity ration parameter. From these figures, fluid speed is enhanced via more stretching of melting surface. The large stretching of the surface is the reason for the more enhancements in fluid motion. Hence, *δ* is favorable to achieve the enhancement in motion of the fluid particles.

**Assessments of heat energy via physical parameters:** the characterization of the thermal energy mechanism against the variation in *Pr*, *α*1, Y<sup>1</sup> and *R* is conducted in Figures 7–10. In Figure 7, the role of second grade fluid number is considered. The mechanism of thermal energy is reduced using higher values of *R*. The thickness of thermal boundary is also decreased. The characterization of heat energy against the distribution of very small number (Y1) based on thermal conductivity is measured by Figure 8. The production of heat energy becomes higher via higher values of Y1. Hence, Y<sup>1</sup> plays an essential role in the enhancement of maximum production of energy. The performance of Prandtl number versus thermal energy is visualized by Figure 9. The temperature profile is noticed as reducing the role versus the variation of *Pr*. This reducing role happens because of the definition of the Prandtl number. According to the definition of *Pr*, MBL is enhanced while the thickness of TBL is increased. Meanwhile, the production of heat energy is reduced with respect to large values of *Pr*. Figure 10 illustrates the trend of *α*<sup>1</sup> versus the temperature profile and *α*<sup>1</sup> is denoted as a relation time parameter. It is investigated how *α*<sup>1</sup> has appeared due to Cattaneo heat flux theory in the energy equation. The thermal energy is reduced due to attendance of non-Fourier's concept whereas thermal energy for the case of Fourier's law is higher than thermal energy for the case of the non-Fourier's concept.

**Figure 2.** Behavior of *f* 0 (*η*) for *M* when *R* = 0.1, *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *δ* = 0.8.

**Figure 3.** Behavior of *g* 0 (*η*) for *M* when *R* = 0.1, *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *δ* = 0.8.

**Figure 4.** Behavior of *f* 0 (*η*) for *R* when *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *M* = 0.1, *δ* = 0.8.

**Figure 5.** Behavior of *g* 0 (*η*) for *R* when *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *M* = 0.1, *δ* = 0.8.

**Figure 6.** Behavior of *g* 0 (*η*) for *δ* when *R* = 0.1, *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *M* = 0.1.

**Figure 7.** Behavior of *θ*(*η*) for *R* when *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>1</sup> = 0.3 = *γ*2, *M* = 0.1, *δ* = 0.8.

**Figure 8.** Behavior of *θ*(*η*) for *γ*<sup>1</sup> when *R* = 0.1, *Sc* = 0.6, *Pr* = 1.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Figure 9.** Behavior of *θ*(*η*) for *Pr* when *R* = 0.1, *Sc* = 0.6, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Figure 10.** Behavior of *θ*(*η*) for *α*<sup>1</sup> when *R* = 0.1, *Sc* = 0.6, *α*<sup>2</sup> = 0.2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Assessments of mass diffusion via physical parameters:** the distribution in transport of mass diffusion is measured against the variation in *Sc*, *α*<sup>2</sup> and *γ*<sup>2</sup> via Figures 11–14. The mechanism related to mass diffusion versus *Sc* called the Schmidt number is conducted by Figure 11. From this figure, declination is measured in view of the transport of mass using large values of the Schmidt number. The transport of mass becomes slow due to the concept of *Sc*. According to this definition, the diffusion of mass (ratio between momentum and mass diffusivities) has an inverse relation with respect to *Sc*. Hence, the solute slows down considering enlargement in *Sc*. Due to this depreciation occurs in the concentration field. Moreover, the combined enlargement in the Schmidt number and solutal relaxation time lessen the concentration profile. The decreasing graph of the concentration profile is observed versus the values of *Sc*. This decreasing trend and local minimum in solute particles is occurred due to reduction of mass diffusivity. Physical, large values of *Sc* reduce mass diffusivity and this reduction in mass diffusivity is a reason for the local minimum in solute particles. The better performance of solute is the investigated against variation in *R* considering by Figure 12. The range of *R* is 0.1 ≤ *R* ≥ 0.8 is used to obtain maximum production of solute. Figure 13 simulates the variation in mass transport with respect to values of *α*2. The decrement in solute is verified considering large values of *α*2. This physical parameter has the ability to restore a state of equilibrium resulting in the solute becoming slow. Figure 14 captures the role of *γ*<sup>2</sup> versus the diffusion of mass. In this figure, diffusion of mass becomes fast using large values of *γ*2.

**Assessments of Nusselt number, divergent velocity and Sherwood number via physical parameters:** the characterization of surface force, temperature gradient and rate of solute is analyzed considering large values of second grade fluid, stretching ration and time relaxation numbers. These simulations are captured by Tables 2 and 3. From Table 2, it can be noted that an enhancement is addressed in the case of positive values of second grade fluid and time relaxation numbers. Table 2 presents the comparative analysis of current model. The remarkable simulations are verified with published work [7]. Table 4 illustrates the impact of time relaxation numbers against the diffusion of mass. The reparable increasing role of the Sherwood number is measured against positive values of the time relaxation number.

**Figure 11.** Behavior of *φ*(*η*) for *Sc* when *R* = 0.1, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Figure 12.** Behavior of *φ*(*η*) for *R* when *Sc* = 0.6, *α*<sup>1</sup> = 0.2 = *α*2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Figure 13.** Behavior of *φ*(*η*) for *α*<sup>2</sup> when *R* = 0.1, *Sc* = 0.6, *α*<sup>1</sup> = 0.2, *γ*<sup>2</sup> = 0.3, *M* = 0.1, *δ* = 0.8.

**Figure 14.** Behavior of *φ*(*η*) for *γ*<sup>2</sup> when *R* = 0.1, *Sc* = 0.6, *α*<sup>1</sup> = 0.2 = *α*2, *M* = 0.1, *δ* = 0.8.



**Table 3.** Comparative investigation for heat transfer rate against *α*<sup>1</sup> .


**Table 4.** Comparative investigation for mass transfer rate against *α*2.


#### **5. Conclusions and Key Findings of the Investigation Performed**

The physical occurrence of solute, heat energy and flow phenomena in a second grade liquid were visualized passing a 3D melting moving surface. The theory of nonFourier's law was imposed in the current flow model inserting variable properties. The current complex model was simulated with the help of an analytical scheme. The main consequences of the current problem are addressed below:


**Author Contributions:** Conceptualization, M.S. and U.N.; methodology, M.S.; software, U.N.; validation, R.A.E.-N., M.M.S. and O.B.; formal analysis, O.B.; investigation, P.T.; resources, H.A.; data curation, U.N.; writing—original draft preparation, M.S.; writing—review and editing, M.S.; visualization, U.N.; supervision, M.S.; project administration, M.S.; funding acquisition, R.A.E.-N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data used to support this study are included in the manuscript.

**Conflicts of Interest:** The authors declare that they have no conflicts of interest.


## **References**

