Fractal Numerical Investigation of Mixed Convective Prandtl-Eyring Nanofluid Flow with Space and Temperature-Dependent Heat Source

Abstract

An explicit computational scheme is proposed for solving fractal time-dependent partial differential equations (PDEs). The scheme is a three-stage scheme constructed using the fractal Taylor series. The fractal time order of the scheme is three. The scheme also ensures stability. The approach is utilized to model the time-varying boundary layer flow of a non-Newtonian fluid over both stationary and oscillating surfaces, taking into account the influence of heat generation that depends on both space and temperature. The continuity equation of the considered incompressible fluid is discretized by first-order backward difference formulas, whereas the dimensionless Navier–Stokes equation, energy, and equation for nanoparticle volume fraction are discretized by the proposed scheme in fractal time. The effect of different parameters involved in the velocity, temperature, and nanoparticle volume fraction are displayed graphically. The velocity profile rises as the parameter I grows. We primarily apply this computational approach to analyze a non-Newtonian fluid’s fractal time-dependent boundary layer flow over flat and oscillatory sheets. Considering spatial and temperature-dependent heat generation is a crucial factor that introduces additional complexity to the analysis. The continuity equation for the incompressible fluid is discretized using first-order backward difference formulas. On the other hand, the dimensionless Navier–Stokes equation, energy equation, and the equation governing nanoparticle volume fraction are discretized using the proposed fractal time-dependent scheme.

1. Introduction

Nanofluids are potential candidates for heat transfer and flow and have been studied profoundly in recent years. Nanofluids are applied in various industries, e.g., in cooling systems, energy storage, and catalysis. Nanofluids are fluids containing nanoparticles ranging from 1 to 100 nm, so they are engineered fluids. This is a mixture of nanoparticles of a specific size tuning with a base fluid.

Nanoparticles also enhance some properties of the fluid (effective property), such as the thermal conductivity of the fluid. The Prandtl-Eyring model is one of the most used models for modelling the flow of these nanofluids.

The nanoparticle in the nanofluid acts somewhat like a gas. The behavior of the nanoparticles is considered in this model. The Prandtl-Eyring model is used to model the flow of the nanofluids and predict their heat transfer performance in various situations, geometries, or configurations.

The Prandtl-Eyring model considers the phenomenon of slipping among the fluid and nanoparticles. But it is not without limitations; first, the fluid-time thermal conductivity is considered constant. Secondly, the behavior and the effects of the size and shape of the nanoparticles on the properties of the fluid are not taken into account. Thus, to overcome these limitations and assumptions, the generalization of the Prandtl-Eyring model can give rise to as many new models, such as the tempered fractional Prandtl-Eyring model, tempered Prandtl-Eyring model, or fractal Prandtl-Eyring model.

When an aircraft wing or a pipe has a liquid layer in contact with its surface, and when the liquid sheet is thin, the limit of fluid mechanics comes into action. The liquid can separate if the shear force affects it in the boundary layer.

This work aims to improve the complex interplay between many components that govern nanoscale fluid flow, heat transmission, and particle dispersion. The primary objective of this undertaking is to gain a comprehensive understanding of and improve the complex interplay among several factors that influence fluid dynamics, heat transmission, and particle dispersion at the nanoscale. When the fluid stays in contact with the surface, the individual layer speeds up the fluid from the optimum limit to zero at the boundary layer. For example, the limits on the trailing edge of the airplane wing are tiny and denser. Along the front or at the upstream of these, the boundaries of the limits thicken. The concept of boundary layers was suggested by Prandtl in 1904 to explain the flow behavior of the viscous fluid near a resolved barrier [1]. Using the finite Reynolds number flow, Prandtl made and deduced the boundary layer equations by using the Navier Stoke equations. As requested, the boundary layer theory equations are the compulsory and necessary simplification of the original Navier–Stokes derivative.

To develop equations on this phenomenon, the vast region over and around the boundary layer and a smaller area where the boundary layer flows are studied closely and built. Various product-specific boundary conditions about the physical model are used to solve the boundary layer equations.

The heat transport in the thermal boundary layer can be reduced by using magnetic and Peclet numbers in the case of magnetohydrodynamic (MHD) flow of fluid with gyrotactic microorganisms [2]. The magnified thermal boundary layer thickness is noticed when a shrunken fluid flows to the growing porous wedge of the Casson liquid as the convective heat transfers [3]. A few other experiments include the influences of various physical parameters on the boundary layer flow [4,5,6].

Many systems use nanofluids, defined as nanomaterials dispersed in a pure fluid, because of their improved specification [7]. Nanofluids have an impressive advantage over more traditional fluids in terms of enhanced thermal conductivity [8]. The interaction between nanoparticles and base fluid, as well as other factors, dramatically affects the viscosity of nanofluids [9]. Although the non-Newtonian behavior of nanofluids has been widely documented, different authors have noted that they behave like Newtonian fluids [10]. Coating wire and blades, processing molten plastic, dyeing textiles, moving biological fluids, processing foods and slurries, and certain petroleum fluids are all examples of practical applications of non-Newtonian behaviors. Power law, micropolar, Reiner-Philippoff, viscoelastic, Casson, Carreau, Giesekus, Prandtl, Prandtl-Eyring, and Powell-Eyring models are among those that can be anticipated in this respect [11]. To achieve the required behavior, these models prioritize the momentum conservation equations to a greater extent. Each model provides a mathematical representation of the relationship between shear stress and deformation rate. The function of sine hyperbolic is Prandtl, while Prandtl-Eyring is a sine inverse function [12,13]. The power law model describes the relationship as nonlinear [14]. A micropolar Prandtl fluid in a porous stretched sheet scenario was investigated by Sajid et al. [15]. We assumed a chemical reaction was occurring in the medium and that the heat source was temperature-dependent. A numerical investigation of power law nanofluid flow on a porous plate was carried out by Maleki et al. [16]. Despite their expectations, they found that using Newtonian nanofluid did not enhance heat transmission.

On the other hand, increasing heat transport was a crucial function of the non-Newtonian one. The phenomenon of Prandtl-Eyring fluid movement through a sensor surface under conditions of magnetic force was studied by Shankar and Naduvinamani [17]. An increase in the velocity field and a decrease in the temperature profile were seen as a result of adding magnetic parameters. Heat transport decreased as the Prandtl number increased simultaneously. Al-Kaabi and Al-Khafajy [18] studied the effects of temperature and concentration changes on the Prandtl-Eyring fluid heat transfer in a porous medium. The experiment was performed by Hayat et al. [19]. According to them, Prandtl-Eyring liquid has been subjected to an extending piece with gyratory microorganisms. It has been proved that by increasing the melting parameters, the velocity increases and decreases the temperature.

Therefore, extending a surface is a well-known and common exercise in many industries, such as extrusion, cooling metallic plates, glass-blown products, ceramics production, rubbers pouring in fiberglass, and many more. One theory that sheds light on the underlying scientific facts is boundary layer flow and heat transfer [20]. Practical difficulties encountered in engineering applications often involve nonlinear equations. Using the Keller-box method [21] can cure such problems. The Prandtl-Eyring liquid flow on a strained sheet was prescribed by Munjam et al. [22]. Munjam et al. [22] also discovered a new method to solve these problems. They then separated their findings from those of the Keller-box system. Their results showed that the velocity also rises with increased liquid parameters. The analysis showed that the velocity increases with a rise in the fluid parameter. They also discovered that, in contrast to the viscous fluid, the Prandtl-Eyring fluid causes a higher gross velocity value.

In their study on the entropy generation of Casson nanofluid, Jamshed et al. [23] utilized the Tiwari and Das model in conjunction with the Keller box approach to solving ordinary differential equations. After incorporating $Cu$ and $Ti{O}_2$ nanoparticles into two different methanol-based nanofluids, the performance of the $Cu$ nanofluids was shown to be superior. Similar concepts and methods were employed in [24] for the base fluid of engine oils and nanoparticles. Their research led them to conclude that increasing the Reynolds and Brinkman numbers would improve entropy generation. Shear rate enhancement was also observed as the concentration of nanofluids increased. Abdelmalek et al. [25] have stated that a Prandtl-Eyring nanofluid can impact Brownian motion and thermophoretic force on a stretched surface. The negative effect of magnetic force on momentum has been demonstrated. On the other hand, the thermal energy was increased via thermophoretic force and Brownian motion.

The growth of various industries, especially in terms of thermal properties, has generated interest in many forms of nanofluid. The applications of nanofluid flow in different sectors are as follows:

• Cooling systems: Nanofluids have demonstrated potential in improving the thermal efficiency of cooling systems, specifically in electronic devices and microchips. Efforts are currently being made to commercialize nanofluids in order to enhance heat removal efficiency in these applications.

• Energy storage: It is worth noting that nanofluids have the potential to enhance the thermal efficiency of energy storage systems, including fuel cells. Efforts are being made to utilize the improved thermal conductivity of nanofluids to optimize energy storage and conversion processes.

• Catalysis: The unique features of nanofluids make them highly promising for catalytic applications, where they can significantly improve reaction speeds and selectivity. Pharmaceutical, chemical synthesis, and environmental remediation are some of the potential applications of nanofluids that are the focus of current research.

• Heat exchangers. Because of their exceptional thermal qualities, nanofluids are finding more and more uses as coolants in heat exchangers, where they greatly enhance performance and efficiency. Improvements in the preparation and use of nanofluids for different heat exchange purposes are the subject of continuing research.

• Vehicle thermal management: Vehicles, such as hybrid electric vehicles and chilled engines, can potentially benefit from nanofluids for heat management. Improving vehicle cooling and thermal management systems is the focus of research efforts. Nanofluids are being considered as potential solutions.

In general, nanofluids are seen as the next big thing in heat transfer fluids since they outperform conventional fluids in terms of thermal conductivity and convective heat transfer coefficient.

The use of spatial fractional derivatives in Fourier’s law modelling has recently been explored in several field studies to improve the precision of physical problem simulations [26,27]. When trying to simulate the behavior of nanofluids, the first-order derivative of the temperature, which is a part of classical Fourier’s law that describes the process of heat conduction, is insufficiently precise. According to recent studies, nanofluid flow modelling can be more accurate using Fourier’s law with spatial-fractional derivatives [28,29]. A fractional order gradient, expressed as the Fourier’s law heat flux relative to the nanofluid, is computed in his model using $q\propto {\nabla }^{\mathsf{\lambda }}T(X,Y)$ where $\lambda$ is a non-integer order.

Caputo and the Riemann-Liouville (R-L) derivatives are popular fractional derivatives [30,31]. Asjad et al. recently applied time-fractional derivatives to nanofluid convection flow between parallel plates [32] and maxwell fluid heat transfer over a vertical surface [33,34,35]. The boundary conditions impact the outcomes more than the initial conditions due to the nonlinear nature of the space derivative terms (e.g., momentum and energy) in the transport equations. Recent research on boundary layer problems has explored fractional derivatives from various angles, as shown in Table 1.

By adding two additional variables to the original equation, which converts it into a system of low-order equations, a completely discrete difference scheme is obtained [43] for a diffusion-wave system. Their innovative approach to tackling time-fractional initial-boundary value problems is detailed in [44]. These problems involve summing Caputo derivatives with orders ranging from 0 to 1. A new linearized transformed $L_1$ Galerkin finite element approach is presented for numerical solutions of the multi-dimensional time fractional Schrödinger equations [45]. The fully discrete scheme’s conditionally optimum error estimates are demonstrated. In reference [46], the fractional boundary conditions multi-term time-space diffusion equation is examined. The standard and shifted Grundwald Letnikov formula is used to approximate the fractional derivative in space, whereas the $L_2-1_{\sigma }$ formula is used to estimate the fractional derivative in time.

We have proposed an efficient three-stage scheme for the time fractal Prandtl-Eyring nanofluid flow with space and temperature-dependent heat sources. The Prandtl-Eyring model is limited in its ability to account for the influence of nanoparticle temperature, shape, and size on the characteristics of the fluid. Utilizing computational techniques, significant insights into the system were derived from investigations into the complex dynamics of a non-Newtonian fluid’s fractal time-dependent boundary layer flow over both flat and oscillatory sheets. By discretizing important equations controlling fluid behavior, such as the continuity equation, the dimensionless Navier–Stokes equation, and the energy and nanoparticle volume fraction equations, the proposed scheme demonstrates its flexibility and resilience in dealing with complex fluid mechanics situations.

Reasons for choosing the method: A fractal Runge–Kutta scheme is proposed. The scheme is explicit and provides the third accurate solution in fractal time. The technique can also be used to find the solutions of classical partial differential equations. It is among the most efficient numerical schemes among other third-order schemes. It does not require linearization, so it can find solutions to nonlinear differential equations without using any other iterative method. Since there is not much work on fractal numerical schemes in the literature, a fractal scheme is proposed. However, the scheme can be used only with time fractal derivatives, in which a classical rise in the dependent variable is used, and fractal variation is considered with the run. So, the considered fractal derivative is just a classical or standard rise over a fractal run.

Pros and cons of the method: The pros and cons can be seen from the constructed Table 2 comparisons of the three schemes. The proposed scheme produces less error than the other two but has a low numerical convergence order. Still, theoretically, it gives a solution of third-order accuracy, whereas the other two methods give the solution of first and second-order accuracy. Nevertheless, the approach may provide difficulties concerning its implementation, substantial computing expenses, restricted validation using empirical data, and susceptibility to parameter selections.

Here is a concise summary of the paper. Two definitions are offered, followed by the proposal of a fractal scheme. The subsequent section presents an analysis of the stability of the scheme for a scalar fractal time partial differential equation, as well as an examination of convergence for a system of time fractal convection diffusion issues. The scheme is also utilized in a fluid dynamics problem, where the results are analyzed and shown in various graphs. Additionally, a comparison of the scheme is conducted.

Definition 1 

([47]). The fractal derivative of a function $g\left(t\right)$ with respect to fractal measure, $t$ is defined as

$\frac{\partial f}{\partial t^{\alpha }}=\underset{t_1\to t}{lim}\frac{g\left(t_1\right)-g\left(t\right)}{t_1^{\alpha }-t^{\alpha }},\alpha >0.$

Definition 2 

([48]). The fractal Taylor series expansion for the function $g\left(t\right)$ is given as

$g\left(t\right)=\sum _{m=0}^{\mathrm{\infty }}\frac{{\left(D^{\alpha }\right)}^{m}g\left(t\right)}{m!},$

where ${\left(D^{\alpha }\right)}^m$, $mth$ order fractal derivative of $g\left(t\right)$ and $0<\alpha <1$.

2. Proposed Fractal Numerical Scheme

This work proposes a fractal computational scheme to solve fractal partial differential equations. The scheme is a three-stage scheme that utilizes the two predictors and one corrector stage. It only discretized the time-dependent terms in a given fractal partial differential equation. This scheme divides the whole-time length into different time levels. It is noted that equal time stepping is used for each stage. However, time steps can be different in all three stages, i.e., stage one utilizes a one-time step, stage two uses another time step, and the third stage employs a time step different from the first two time steps. To apply the scheme, consider the fractal partial differential equation of the following form:

(1)$\frac{\partial v}{\partial t^{\alpha }}=G\left(v,\frac{\partial v}{\partial x},\frac{\partial v}{\partial y},\frac{{\partial }^{2}v}{\partial y^2}\right),$

where $0<\alpha \le 1$. Subject to initial and boundary conditions given as follows:

(2)$\left.\begin{array}{c}v\left(x,y,0\right)=0\\ v\left(0,y,t\right)={\beta }_{1}v\left(x,0,t\right)={\beta }_2\\ v\left(x,L,t\right)={\beta }_3\end{array}\right\},$

where ${\beta }_{i}i=\mathrm{1,2},3$ are constant.

The first stage of the proposed scheme is described as follows:

(3)${\overline{v}}_{i,j}^{n+1}=v_{i,j}^n+\mathsf{\Delta }{t}_1{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n,$

where $\mathsf{\Delta }{t}_1=\frac{1}{3}(t_{n+1}^{\alpha }-t_n^{\alpha })$ is a time step.

The second stage of the scheme is given as follows:

(4)${\stackrel{̿}{v}}_{i,j}^{n+1}=\frac{1}{3}\left(2v_{i,j}^n+{\overline{v}}_{i,j}^{n+1}\right)+\mathsf{\Delta }{t}_2{\left.\frac{\partial \overline{v}}{\partial t^{\alpha }}\right|}_{i,j}^{n+1}=v_{i,j}^n+\left(\frac{1}{3}\mathsf{\Delta }{t}_1+\mathsf{\Delta }{t}_2\right){\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n$

Stages (3) and (4) are predictor stages. The corrector stage can be expressed as follows:

(5)$v_{i,j}^{n+1}=a{v}_{i,j}^n+b{\overline{v}}_{i,j}^{n+1}+c{\stackrel{̿}{v}}_{i,j}^{n+1}+d\mathsf{\Delta }t{\left.\frac{\partial \stackrel{̿}{v}}{\partial t^{\alpha }}\right|}_{i,j}^{n+1},$

where $a,b,c,$ and $d$ are parameters to be determined later.

By using the first and second stages (3) and (4) in the third stage (5), it results in the following:

(6)$\begin{array}{ll}{v}_{i,j}^{n+1}& =a{v}_{i,j}^n+b{v}_{i,j}^n+b\mathsf{\Delta }{t}_1{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+c{v}_{i,j}^n+c\left(\frac{1}{3}\mathsf{\Delta }{t}_1+\mathsf{\Delta }{t}_2\right){\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+c\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2{\left(\left.\frac{\partial v}{\partial t^{\alpha }}\right|\right)}^2{v}_{i,j}^n\\ & +d\mathsf{\Delta }t\left\{{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+\left(\frac{1}{3}\mathsf{\Delta }{t}_1+\mathsf{\Delta }{t}_2\right){\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n+\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^3{v}_{i,j}^n\right\},\end{array}$

Consider the fractal Taylor series expansion $v_{i,j}^{n+1}$ as follows:

(7)$v_{i,j}^{n+1}=v_{i,j}^n+\mathsf{\Delta }t\left(\frac{\partial }{\partial t^{\alpha }}\right)v_{i,j}^n+\frac{{\left(\mathsf{\Delta }t\right)}^2}{2}{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n+\frac{{\left(\mathsf{\Delta }t\right)}^3}{6}{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^3{v}_{i,j}^n+\dots$

Substituting fractal Taylor series expansion (7) into Equation (6) results in the following:

(8)$\begin{array}{ll}{v}_{i,j}^n+\mathsf{\Delta }t{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n& +\frac{{\left(\mathsf{\Delta }t\right)}^2}{2}{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n+\frac{{\left(\mathsf{\Delta }t\right)}^3}{6}{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^3{v}_{i,j}^n+O\left({\left(\mathsf{\Delta }t\right)}^4\right)\\ & =a{v}_{i,j}^n+b{v}_{i,j}^n+b\mathsf{\Delta }{t}_1{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+c{v}_{i,j}^n+c\left(\frac{1}{3}\mathsf{\Delta }{t}_1+\mathsf{\Delta }{t}_2\right){\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+c\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n\\ & +d\mathsf{\Delta }t\left\{{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n+\left(\frac{1}{3}\mathsf{\Delta }{t}_1+\mathsf{\Delta }{t}_2\right){\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n+\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^3{v}_{i,j}^n\right\},\end{array}$

By equating coefficients of $v_{i,j}^n,{\left.\frac{\partial v}{\partial t^{\alpha }}\right|}_{i,j}^n,{\left(\frac{\partial }{\partial t^{\alpha }}\right)}^2{v}_{i,j}^n$ and ${\left(\frac{\partial }{\partial t^{\alpha }}\right)}^3{v}_{i,j}^n$ on both sides of Equation (8) it yields the following:

(9)$\left.\begin{array}{c}a+b+c=1\\ b\Delta t_1+\frac{c}{3}\Delta t_1+c\Delta t_2+d\Delta t=\Delta t\\ c\Delta t_1\Delta t_2+\frac{d}{3}\Delta t\Delta t_1+d\Delta t\Delta t_2=\frac{{\left(\mathsf{\Delta }t\right)}^2}{2}\\ d\Delta t\Delta t_1\Delta t_2=\frac{{\left(\mathsf{\Delta }t\right)}^3}{6}\end{array}\right\},$

By solving the system of Equation (9), the values of unknown can be written as follows:

(10)$\left.\begin{array}{c}a=\frac{-2\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_1^2-12\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2+18\mathsf{\Delta }{t}^2\mathsf{\Delta }{t}_1^2\mathsf{\Delta }{t}_2+9\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_2^2-27\mathsf{\Delta }{t}^2\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2^3+54\mathsf{\Delta }t\mathsf{\Delta }{t}_1^2\mathsf{\Delta }{t}_2^2-54\mathsf{\Delta }{t}_1^3\mathsf{\Delta }{t}_2^2}{54\mathsf{\Delta }{t}_1^3\mathsf{\Delta }{t}_2^2},\\ b=\frac{\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_1^2+3\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2+9\mathsf{\Delta }{t}^2\mathsf{\Delta }{t}_1^2\mathsf{\Delta }{t}_2-9\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_2^2+27\mathsf{\Delta }{t}^2\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2^2-54\mathsf{\Delta }t\mathsf{\Delta }{t}_1^2\mathsf{\Delta }{t}_2^2}{54\mathsf{\Delta }{t}_1^3\mathsf{\Delta }{t}_2^2},\\ c=\frac{-\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_1+32\mathsf{\Delta }{t}^3\mathsf{\Delta }{t}_2-9\mathsf{\Delta }{t}^2\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2}{18\mathsf{\Delta }{t}_1^2\mathsf{\Delta }{t}_2^2},\\ d=\frac{\mathsf{\Delta }{t}^2}{6\mathsf{\Delta }{t}_1\mathsf{\Delta }{t}_2}.\end{array}\right\}$

Therefore, the proposed scheme for fractal discretization of Equation (1) with $G={\alpha }_1\frac{\partial v}{\partial x}+{\alpha }_2\frac{\partial v}{\partial y}+{\alpha }_3\frac{{\partial }^{2}v}{\partial y^2}$ and second-order space discretization is written as follows:

(11)${\overline{v}}_{i,j}^{n+1}=v_{i,j}^n+\mathsf{\Delta }{t}_1\left({\alpha }_1{\delta }_x{v}_{i,j}^n+{\alpha }_2{\delta }_y{v}_{i,j}^n+{\alpha }_3{\delta }_y^2{v}_{i,j}^n\right).$

(12)${\stackrel{̿}{v}}_{i,j}^{n+1}=\frac{1}{3}\left(2v_{i,j}^n+{\overline{v}}_{i,j}^{n+1}\right)+\mathsf{\Delta }{t}_2\left({\alpha }_1{\delta }_x{\overline{v}}_{i,j}^{n+1}+{\alpha }_2{\delta }_y{\delta }_x{\overline{v}}_{i,j}^{n+1}+{\alpha }_3{\delta }_y^2{\overline{v}}_{i,j}^{n+1}\right).$

(13)$v_{i,j}^{n+1}=a{v}_{i,j}^n+b{\overline{v}}_{i,j}^{n+1}+c{\stackrel{̿}{v}}_{i,j}^{n+1}+d\mathsf{\Delta }t\left({\alpha }_1{\delta }_x{\stackrel{̿}{v}}_{i,j}^{n+1}+{\alpha }_2{\delta }_y{\stackrel{̿}{v}}_{i,j}^{n+1}+{\alpha }_3{\delta }_y^2{\stackrel{̿}{v}}_{i,j}^{n+1}\right).$

where ${\delta }_x{v}_{i,j}^n=\frac{v_{i+1,j}^n-v_{i-1,j}^n}{2\mathsf{\Delta }x},{\delta }_y{v}_{i,j}^n=\frac{v_{i,j+1}^n-v_{i,j-1}^n}{2\mathsf{\Delta }y},{\delta }_y^2{v}_{i,j}^n=\frac{v_{i,j+1}^n-2v_{i,j}^n+v_{i,j-1}^n}{{\left(\mathsf{\Delta }y\right)}^2}$.

3. Stability Analysis

The Fourier series analysis can be used for finite difference methods to find the stability conditions of linear differential equations. The technique transforms the difference equation into trigonometric equations. The analysis is utilized for linear partial differential equations, while it determines the precise stability criteria for nonlinear partial differential equations. Because for nonlinear partial differential equations, it is used to linearize PDEs. To apply this analysis, consider the following transformations:

(14)$\left.\begin{array}{c}{v}_{i,j}^n=E^n{e}^{iJ{\psi }_1}{e}^{jI{\psi }_2}\\ {\overline{v}}_{i,j}^{n+1}={\overline{E}}^{n+1}{e}^{iI{\psi }_1}{e}^{jI{\psi }_2}\\ {\stackrel{̿}{v}}_{i,j}^{n+1}={\stackrel{̿}{E}}^{n+1}{e}^{iI{\psi }_1}{e}^{jI{\psi }_2}\end{array}\right\},$

where $I=\sqrt{-1}$ and ${\psi }_1=\overline{p}∆x$ and ${\psi }_2=\overline{q}∆y$ are two-phase angles and $\overline{p}$ and $\overline{q}$ are two wave numbers, respectively, in $x$ and $y$ directions, i.e.$,$ ${\mathsf{\lambda }}_{\mathrm{x}}=\frac{2\mathsf{\pi }}{\overline{p}}$ and ${\mathsf{\lambda }}_{\mathrm{y}}=\frac{2\mathsf{\pi }}{\overline{q}}$ are wavelengths, and $E^n$ is the amplitude at time level $n$. By substituting the relevant transformations from Equation (14) into the first step of the proposed scheme Equation (11), it produces a result.

(15)$\begin{array}{ll}{\overline{E}}^{n+1}{e}^{iI{\psi }_1}{e}^{jI{\psi }_2}& =E^n{e}^{iI{\psi }_1}{e}^{jI{\psi }_2}\\ & +\mathsf{\Delta }{t}_1\{{\alpha }_1\left(\frac{e^{\left(i+1\right)I{\psi }_1}{e}^{jI{\psi }_2}-e^{\left(i-1\right)I{\psi }_1}{e}^{jI{\psi }_2}}{2\mathsf{\Delta }x}\right)E^n\\ & +{\alpha }_2\left(\frac{e^{iI{\psi }_1}{e}^{\left(j+1\right)I{\psi }_2}-e^{iI{\psi }_1}{e}^{\left(j-1\right)I{\psi }_2}}{2\mathsf{\Delta }y}\right)E^n\\ & +{\alpha }_3\left(\frac{e^{iI{\psi }_1}{e}^{(j+1)I{\psi }_2}-2e^{iI{\psi }_1}{e}^{jI{\psi }_2}+e^{iI{\psi }_1}{e}^{\left(j-1\right)I{\psi }_2}}{{\left(\mathsf{\Delta }y\right)}^2}\right)E^n\},\end{array}$

Simplifying Equation (15) yields

${\overline{E}}^{n+1}=E^n+\mathsf{\Delta }{t}_1\left\{{\alpha }_1\left(\frac{e^{I{\psi }_1}-e^{-I{\psi }_1}}{2\mathsf{\Delta }x}\right)E^n+{\alpha }_2\left(\frac{e^{I{\psi }_2}-e^{-I{\psi }_2}}{2\mathsf{\Delta }y}\right)E^n+{\alpha }_3\left(\frac{e^{I{\psi }_2}-2+e^{-I{\psi }_2}}{{\left(\mathsf{\Delta }y\right)}^2}\right)E^n\right\}$

(16)${\overline{E}}^{n+1}=E^n+\mathsf{\Delta }{t}_1\left\{{\alpha }_1\frac{Isin{\psi }_1}{\mathsf{\Delta }x}+{\alpha }_2\frac{Isin{\psi }_2}{\mathsf{\Delta }y}+{\alpha }_3\frac{(2cos{\psi }_2-2)}{{\left(\mathsf{\Delta }y\right)}^2}\right\}{E}^n.$

Let ${\overline{c}}_1=\frac{{\alpha }_1\mathsf{\Delta }{t}_1}{\mathsf{\Delta }x},{\overline{c}}_2=\frac{{\alpha }_2\mathsf{\Delta }{t}_1}{\mathsf{\Delta }y},{\overline{d}}_1=\frac{{2\alpha }_3\mathsf{\Delta }{t}_1}{{\left(\mathsf{\Delta }y\right)}^2},$ then Equation (16) can be written as follows:

(17)${\overline{E}}^{n+1}=E^n+\left\{{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2+{\overline{d}}_1(cos{\psi }_2-1)\right\}{E}^n.$

Similarly, the second stage can be transformed as follows:

(18)$\begin{array}{ll}{\overline{\overline{E}}}^{n+1}& =\frac{1}{3}\left(2E^n+{\overline{E}}^{n+1}\right)+\left\{{\widehat{c}}_{1}Isin{\psi }_1+{\widehat{c}}_{2}Isin{\psi }_2+{\widehat{d}}_1(cos{\psi }_2-1)\right\}{\overline{E}}^{n+1}\\ & =\frac{2}{3}{E}^n+\left(\frac{1}{3}+{\widehat{c}}_{1}Isin{\psi }_1+{\widehat{c}}_{2}Isin{\psi }_2+{\widehat{d}}_1\left(cos{\psi }_2-1\right)\right){\overline{E}}^{n+1}\\ & =\frac{2}{3}{E}^n\\ & +\left(\frac{1}{3}+{\widehat{c}}_{1}Isin{\psi }_1+{\widehat{c}}_{2}Isin{\psi }_2+{\widehat{d}}_1\left(cos{\psi }_2-1\right)\right)(1+{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2+{\overline{d}}_1(cos{\psi }_2\\ & -1))E^n.\end{array}$

where ${\widehat{c}}_1={\alpha }_1\frac{\mathsf{\Delta }{t}_2}{\mathsf{\Delta }x},{\widehat{c}}_2={\alpha }_2\frac{\mathsf{\Delta }{t}_2}{\mathsf{\Delta }y},{\widehat{d}}_1=2{\alpha }_3\frac{\mathsf{\Delta }{t}_2}{{\left(\mathsf{\Delta }y\right)}^2}$.

The third stage of the scheme can be expressed as follows:

(19)$\begin{array}{ll}{E}^{n+1}={aE}^n& +b\left\{1+{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2+{\overline{d}}_1\left(cos{\psi }_2-1\right)\right\}{E}^n\\ & +\left\{c+c_{1}Isin{\psi }_1+c_{2}Isin{\psi }_2+d_1\left(cos{\psi }_2-1\right)\right\}\{\frac{2}{3}\\ & +\left(\frac{1}{3}+{\widehat{c}}_{1}Isin{\psi }_1+{\widehat{c}}_{2}Isin{\psi }_2+{\widehat{d}}_1\left(cos{\psi }_2-1\right)\right)(1+{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2\\ & +2{\overline{d}}_1\left(cos{\psi }_2-1\right)\left)\right\}{E}^n,\end{array}$

The amplification factor, in this case, is

(20)$\begin{array}{ll}{\left|\frac{E^{n+1}}{E^n}\right|}^2=& |a+b\left\{1+{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2+{\overline{d}}_1\left(cos{\psi }_2-1\right)\right\}\\ & +\left\{c+c_{1}Isin{\psi }_1+c_{2}Isin{\psi }_2+d_1\left(cos{\psi }_2-1\right)\right\}\{\frac{2}{3}\\ & +\left(\frac{1}{3}+{\widehat{c}}_{1}Isin{\psi }_1+{\widehat{c}}_{2}Isin{\psi }_2+{\widehat{d}}_1\left(cos{\psi }_2-1\right)\right)(1+{\overline{c}}_{1}Isin{\psi }_1+{\overline{c}}_{2}Isin{\psi }_2\\ & +2{\overline{d}}_1\left(cos{\psi }_2-1\right)\left)\right\}{|}^2<1.\end{array}$

Therefore, the scheme will be stable if it satisfies the conditions (20).

The stability of the suggested scheme for scalar partial differential equations has been demonstrated. Now, we will provide the convergence of the method for the system of partial differential equations. To do so, consider the following vector matrix equation:

(21)$\frac{\partial \mathit{g}}{\partial t^{\alpha }}=A_1\frac{\partial \mathit{g}}{\partial x}+B_1\frac{\partial \mathit{g}}{\partial y}+C_1\frac{{\partial }^2\mathit{g}}{\partial y^2}+D_1\mathit{g}.$

By applying a proposed scheme on Equation (20), it gives the following:

(22)${\overline{\mathit{g}}}_{i,j}^{n+1}={\mathit{g}}_{i,j}^n+\mathsf{\Delta }{t}_1\left\{A_1{\delta }_x{\mathit{g}}_{i,j}^n+B_1{\delta }_y{\mathit{g}}_{i,j}^n+C\_1{\delta }_y^2{\mathit{g}}_{i,j}^n+D_1{\mathit{g}}_{i,j}^n\right\}.$

(23)${\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}=\frac{1}{3}\left(2{\mathit{g}}_{i,j}^n+{\overline{\mathit{g}}}_{i,j}^{n+1}\right)+\mathsf{\Delta }{t}_2\left(A_1{\delta }_x{\overline{\mathit{g}}}_{i,j}^{n+1}+B_1{\delta }_y{\overline{\mathit{g}}}_{i,j}^{n+1}+C_1{\delta }_y^2{\overline{\mathit{g}}}_{i,j}^{n+1}+D_1{\overline{\mathit{g}}}_{i,j}^{n+1}\right).$

(24)${\mathit{g}}_{i,j}^{n+1}=a{\mathit{g}}_{i,j}^n+b{\overline{\mathit{g}}}_{i,j}^{n+1}+c{\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}+\mathsf{\Delta }t\left(A_1{\delta }_x{\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}+B_1{\delta }_y{\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}+C_1{\delta }_y^2{\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}+D_1{\stackrel{̿}{\mathit{g}}}_{i,j}^{n+1}\right).$

Theorem 1.

The proposed scheme (22)–(24) converges for Equation (21).

Proof. 

To prove this theorem, consider the first stage of the exact scheme.

(25)${\overline{\mathit{G}}}_{i,j}^{n+1}={\mathit{G}}_{i,j}^n+\mathsf{\Delta }{t}_1\left\{A_1{\delta }_x{\mathit{G}}_{i,j}^n+B_1{\delta }_y{\mathit{G}}_{i,j}^n+C_1{\delta }_y^2{\mathit{G}}_{i,j}^n+D_1{\mathit{G}}_{i,j}^n\right\},$

Subtracting Equation (22) from Equation (25), the following is obtained:

(26)${\overline{\mathit{e}}}_{i,j}^{n+1}={\mathit{e}}_{i,j}^n+\mathsf{\Delta }{t}_1\left\{A_1{\delta }_x{\mathit{e}}_{i,j}^n+B_1{\delta }_y{\mathit{e}}_{i,j}^n+C_1{\delta }_y^2{\mathit{e}}_{i,j}^n+D_1{\mathit{e}}_{i,j}^n\right\},$

where ${\overline{\mathit{e}}}_{i,j}^{n+1}={\overline{\mathit{G}}}_{i,j}^{n+1}-{\overline{\mathit{g}}}_{i,j}^{n+1},{\mathit{e}}_{i,j}^n={\mathit{G}}_{i,j}^n-{\mathit{g}}_{i,j}^n$, etc.

Taking the norm on both sides of Equation (25) gives

(27)${\overline{\mathit{e}}}_{i,j}^{n+1}\le {\mathit{e}}_{i,j}^n+\mathsf{\Delta }{t}_1\left\{‖A_1‖\frac{{\mathit{e}}_{i,j}^n}{\mathsf{\Delta }x}+‖B_1‖\frac{{\mathit{e}}_{i,j}^n}{\mathsf{\Delta }y}+‖C_1‖\frac{4{\mathit{e}}_{i,j}^n}{{\left(\mathsf{\Delta }\mathrm{y}\right)}^2}+‖D_1‖{\mathit{e}}_{i,j}^n\right\}.$

Re-write inequality (27) as

(28)${\overline{e}}^{n+1}\le e^n+\left({\overline{P}}_1+{\overline{P}}_2+{\overline{P}}_3+{\overline{P}}_4\right)e^n$

where ${\overline{P}}_1=\frac{‖A_1‖\mathsf{\Delta }{t}_1}{\mathsf{\Delta }x},{\overline{P}}_2=\frac{‖B_1‖\mathsf{\Delta }{t}_1}{\mathsf{\Delta }y},{\overline{P}}_3=\frac{4‖C_1‖}{{\left(\mathsf{\Delta }y\right)}^2}$ and, ${\overline{P}}_4=‖D_1‖$.

Similarly, for the second stage, the following inequality can be obtained.

(29)${\stackrel{̿}{e}}^{n+1}\le \frac{1}{3}\left({2e}^n+{\overline{e}}^{n+1}\right)+\left({\widehat{P}}_1+{\widehat{P}}_2+{\widehat{P}}_3+{\widehat{P}}_4\right){\overline{e}}^{n+1}\le \left\{\frac{2}{3}+\left(\frac{1}{3}+{\widehat{P}}_1+{\widehat{P}}_2+{\widehat{P}}_3+{\widehat{P}}_4\right)\left(1+{\overline{P}}_1+{\overline{P}}_2+{\overline{P}}_3+{\overline{P}}_4\right)\right\}{e}^n.$

where ${\widehat{P}}_1=\frac{\mathsf{\Delta }{t}_2‖A_1‖}{\mathsf{\Delta }x},{\widehat{P}}_2=\frac{\mathsf{\Delta }{t}_2‖B_1‖}{\mathsf{\Delta }y},{\widehat{P}}_3=\frac{\mathsf{\Delta }{t}_2‖C_1‖}{{\left(\mathsf{\Delta }y\right)}^2},{\widehat{P}}_4=\mathsf{\Delta }{t}_2‖D_1‖$.

Similarly, the following inequality can be obtained by subtracting the third stage of the scheme.

(30)$e^{n+1}\le {ae}^n+b{\overline{e}}^{n+1}+c{\stackrel{̿}{e}}^{n+1}+d\left(P_1+P_2+P_3+P_4\right){\stackrel{̿}{e}}^{n+1}+N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right)$

By using inequalities (28) and (29) into inequality (30), which yields

(31)$\begin{array}{ll}{e}^{n+1}\le {ae}^n& +b\left(1+{\overline{P}}_1+{\overline{P}}_2+{\overline{P}}_3+{\overline{P}}_4\right)e^n\\ & +\left\{c+d\left(P_1+P_2+P_3+P_4\right)\right\}\{\frac{2}{3}\\ & +\left(\frac{1}{3}+{\widehat{P}}_1+{\widehat{P}}_2+{\widehat{P}}_3+{\widehat{P}}_4\right)\left(1+{\overline{P}}_1+{\overline{P}}_2+{\overline{P}}_3+{\overline{P}}_4\right)\}{e}^n\\ & +N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right),\end{array}$

Re-write inequality (31) as

(32)$e^{n+1}\le \mu e^n+\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right),$

Let $n=0$ in inequality (32). It gives

(33)$e^1\le N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right),$

Let $n=1$ in inequality (32), it gives

(34)$e^2\le {\mu e}^1+N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right)\le \left(1+\mu \right)N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right),$

If this is continued, then for finite $n$

(35)$e^n\le \left({\mu }^{n-1}+\dots +\mu +1\right)N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right)=\frac{1\left(1-{\mu }^n\right)}{1-\mu }N\left(O\left({\left(\mathsf{\Delta }t\right)}^3,{\left(\mathsf{\Delta }x\right)}^2,{\left(\mathsf{\Delta }y\right)}^2\right)\right).$

By applying the $limitn\to \infty$ in (35), the series $1+\mu +\dots +{\mu }^n+\dots$ becomes an infinite geometric series that will converge if $\left|\mu \right|<1$. □

4. Problem Formulation

Examine the flow of a non-Newtonian fluid over a two-dimensional moving sheet, where the fluid is incompressible, unstable, and laminar. The movement of the sheet generates the flow in the fluid. Let the sheet move with the velocity $u_w.$ The $x^{*}$-axis is placed along the sheet, whereas $y^{*}$-axis is perpendicular to $x^{*}$-axis. The fluid is moving from left to right along $x^{*}$-axis. Temperature and concentration on the sheet are considered higher than ambient temperature and concentration. The governing equations can be written as follows:

(36)$\frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}=0,$

(37)$\frac{\partial u^{*}}{\partial t^{*}}+u^{*}\frac{\partial u^{*}}{\partial x^{*}}+v^{*}\frac{\partial u^{*}}{\partial y^{*}}=\nu \frac{{\mathsf{\Lambda }}_w}{{\overline{C}}_1}\frac{{\partial }^2{u}^{*}}{\partial {y^{*}}^2}+\nu \frac{{\mathsf{\Lambda }}_w}{{2\overline{C}}_1^3}{\left(\frac{\partial u^{*}}{\partial y^{*}}\right)}^2\frac{{\partial }^2{u}^{*}}{\partial {y^{*}}^2}-\frac{\sigma B_{\circ }^2}{\rho }+g\left({\beta }_T\left(T-T_{\infty }\right)+{\beta }_C\left(C-C_{\infty }\right)\right),$

(38)$\frac{\partial T}{\partial t^{*}}+u^{*}\frac{\partial T}{\partial x^{*}}+v^{*}\frac{\partial T}{\partial y^{*}}=\alpha \frac{{\partial }^{2}T}{\partial {y^{*}}^2}+\tau \left(\frac{D_B}{\mathsf{\Delta }C}\frac{\partial T}{\partial y^{*}}\frac{\partial C}{\partial y^{*}}+\frac{D_T}{T_{\infty }}{\left(\frac{\partial T}{\partial y^{*}}\right)}^2\right)+\frac{1}{\rho c_p}{q}^{\prime \prime \prime },$

(39)$\frac{\partial C}{\partial t^{*}}+u^{*}\frac{\partial C}{\partial x^{*}}+v^{*}\frac{\partial C}{\partial y^{*}}=D_B\frac{{\partial }^{2}C}{\partial y^{{*}^2}}+\frac{\mathsf{\Delta }C}{T_{\infty }}{D}_T\frac{{\partial }^{2}T}{\partial {y^{*}}^2}+\frac{D_B}{T_m}{k}_T\frac{{\partial }^{2}T}{\partial {y^{*}}^2}-k_r\left(C-C_{\infty }\right).$

Subject to the initial and boundary conditions

(40)$\left.\begin{array}{c}{u}^{*}=0,v^{*}=0,T=0,C=0when{t}^{*}=0\\ u^{*}=u_w,v^{*}=0,T=T_w,C=C_{w}when{y}^{*}=0\\ u^{*}=0=v^{*}=T=Cwhen{x}^{*}=0\\ u^{*}\to 0,T\to T_{\infty },C\to C_{\infty }when{y}^{*}\to \infty \end{array}\right\},$

where $q^{\prime \prime \prime }=\frac{k{u}_w}{x\nu \rho c_p}\left(A^{*}u\left(T_w-T_{\infty }\right)+B^{*}\left(T-T_{\infty }\right)\right)$ is space- and temperature-dependent internal heat generation. Where ${\mathsf{\Lambda }}_w$ and ${\overline{C}}_1$ are material constants, $g$ is gravity, ${\beta }_T$ represents coefficients of thermal expansion, ${\beta }_C$ represents coefficients of solutal expansion, $\sigma$ is electric conductivity, $B_{\circ }$ is the strength of the magnetic field, $D_B$ represents Brownian motion coefficients, $D_T$ is thermophoresis diffusion, $\mathsf{\Delta }C$ is nanoparticle fraction difference, $T_m$ represents mean fluid temperature, $k_T$ denotes thermal diffusion ratio and $k_r$ is dimensional reaction rate.

By using transformations

(41)$\left.\begin{array}{c}x=\frac{x^{*}}{L},y=\frac{y^{*}}{L},t=\frac{u_w{t}^{*}}{L},u=\frac{u^{*}}{u_w}\\ \theta =\frac{T-T_{\infty }}{T_w-T_{\infty }},\varphi =\frac{C-C_{\infty }}{C_w-C_{\infty }}\end{array}\right\},$

Incorporating transformations (41) into (36)–(40) that yields

(42)$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$

(43)$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{{ϵ}_1}{R_e}\frac{{\partial }^{2}u}{\partial y^2}+{ϵ}_1{ϵ}_2{\left(\frac{\partial u}{\partial y}\right)}^2\frac{{\partial }^{2}u}{\partial y^2}-\frac{H_{\circ }^2}{R_e}u+\frac{{G_{\gamma }}_T}{R_e^2}\theta +\frac{{G_{\gamma }}_C}{R_e^2}\varphi ,$

(44)$\frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\frac{1}{P_r}\frac{1}{R_e}\frac{{\partial }^{2}u}{\partial y^2}+\frac{N_b}{R_e}\frac{\partial \theta }{\partial y}\frac{\partial \varphi }{\partial y}+\frac{N_t}{R_e}{\left(\frac{\partial \theta }{\partial y}\right)}^2+\frac{ϵ}{x}\frac{1}{P_r}\left(A^{*}u+B^{*}\theta \right),$

(45)$\frac{\partial \varphi }{\partial t}+u\frac{\partial \varphi }{\partial x}+v\frac{\partial \varphi }{\partial y}=\frac{1}{S_c}\frac{1}{R_e}\frac{{\partial }^2\varphi }{\partial y^2}+\frac{N_t}{N_b}\frac{1}{S_c{R}_e}\frac{{\partial }^2\theta }{\partial y^2}+\frac{S_r}{R_e}\frac{{\partial }^2\theta }{\partial y^2}-\gamma \varphi .$

Subject to the dimensionless initial and boundary conditions

(46)$\left.\begin{array}{c}u=0,v=0,\theta =0=\varphi ,whent=0\\ u=1,v=0,\theta =1=\varphi ,wheny=0\\ u=0=v,\theta =0=\varphi whenx=0\\ u\to 0,\theta \to 0,\varphi \to 0wheny\to \infty \end{array}\right\},$

where ${ϵ}_1$ represents the parameter I, the dimensionless ratio of material constants indicated as ${ϵ}_1$, is represented by parameter I. The stability and behavior of the fluid flow are significantly influenced by its role. This parameter plays a crucial role in determining the flow characteristics by influencing the strength of fluid motion concerning its viscosity. ${ϵ}_2$ denotes Prandtl-Eyring parameter II, the Prandtl-Eyring parameter, abbreviated as ${ϵ}_2$, is a parameter that represents the dimensionless ratio of the product of the sheet velocity $u_w$, viscosity $\nu$, and characteristic length. The measurement assesses the impact of non-Newtonian phenomena on fluid flow dynamics. The degree of shear-thinning or shear-thickening behavior displayed by the fluid is influenced by Parameter II, which has important consequences for the evolution of the boundary layer and heat transfer. ${G_{\gamma }}_T$ is the thermal Grashof number, ${G_{\gamma }}_C$ is solutal Grashof number, $H_{\circ }$ Hartmann number, $ϵ$ is the dimensionless parameter, $R_e$ is the Reynolds number, $N_t$ represents the thermophoresis parameter, $N_b$ is Brownian motion coefficients, $S_r$ is Soret number, $P_r$ is Prandtl number, $S_c$ is the Schmidt number, and $\gamma$ is the dimensionless parameter, and these are defined as

$$\begin{gathered} {\epsilon }_1=\frac{{\mathsf{\Lambda }}_w}{{\overline{C}}_1},{\epsilon }_2=\frac{u_w\nu }{{2\overline{C}}_1^2{L}^3},{G_{\gamma }}_T=\frac{L^{3}g{\beta }_T\left(T_w-T_{\infty }\right)}{{\nu }^2},{G_{\gamma }}_C=\frac{L^{3}g{\beta }_C\left(C_w-C_{\infty }\right)}{{\nu }^2},H_{\circ }=B_{\circ }L\sqrt{\frac{\sigma }{\rho \nu }},R_e=\frac{L{u}_w}{\nu },N_t \\ =\frac{\tau D_T\left(T_w-T_{\infty }\right)}{\nu T_{\infty }},N_b=\frac{\tau D_B\left(C_w-C_{\infty }\right)}{\nu \mathsf{\Delta }C},S_r=\frac{k_T{D}_B\left(T_w-T_{\infty }\right)}{\nu T_m\left(C_w-C_{\infty }\right)},P_r=\frac{\nu }{\alpha },S_C=\frac{\nu }{D_B},\gamma =\frac{L{k}_r}{u_w} \end{gathered}$$

The skin friction coefficients, local Nusselt number, and local Sherwood number are precisely defined as follows:

(47)$C_f=\frac{{\tau }_w}{\rho u_w^2},$

(48)$N_{u_L}=\frac{L{q}_w}{k\left(T_w-T_{\infty }\right)},$

(49)$S_{h_C}=\frac{L{q}_j}{D}-B\left(C_w-C_{\infty }\right).$

where ${\tau }_w=\mu \frac{{\mathsf{\Lambda }}_w}{{\overline{C}}_1}\frac{\partial u^{\mathrm{*}}}{\partial y^{\mathrm{*}}}-\mu \frac{{\mathsf{\Lambda }}_w}{{6\overline{C}}_1^3}{\left.{\left(\frac{\partial u^{\mathrm{*}}}{\partial y^{\mathrm{*}}}\right)}^3\right|}_{y^{\mathrm{*}}=0},q_w=-k{\left(\frac{\partial T}{\partial y^{\mathrm{*}}}\right)}_{y^{\mathrm{*}}=0},q_j=-D_B{\left(\frac{\partial C}{\partial y^{\mathrm{*}}}\right)}_{y^{\mathrm{*}}=0}$.

By using transformation (41) into (47)–(49), it gives

(50)$C_f={\left(\frac{{\epsilon }_1}{R_e}\frac{\partial u}{\partial y}-\frac{1}{3}{\epsilon }_1{\epsilon }_2{\left(\frac{\partial u}{\partial y}\right)}^3\right)}_{y=0},$

(51)$N_{u_L}=-{\left.\frac{\partial \theta }{\partial y}\right|}_{y=0},$

(52)$S_{h_L}=-{\left.\frac{\partial \varphi }{\partial y}\right|}_{y=0},$

Since fractal derivatives are more generalized than classical ones, they rise over fractal runs. So, it finds the change of the dependent variable over the fractal measure. The time fractal derivative uses only the change of quantities using the fractal shift in time. Physically, fractal change can be considered an irregular interval of time. From a derivation point of view, Navier–Stokes equations can be generalized using fractal derivatives. The proposed scheme can be used for classical and time fractal partial differential equations. The time fractal model is

(53)$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$

(54)$\frac{\partial u}{\partial t^{\alpha }}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{{\epsilon }_1}{R_e}\frac{{\partial }^{2}u}{\partial y^2}+{\epsilon }_1{\epsilon }_2{\left(\frac{\partial u}{\partial y}\right)}^2\frac{{\partial }^{2}u}{\partial y^2}-\frac{H_{\circ }^2}{R_e}u+\frac{{G_{\gamma }}_T}{R_e^2}\theta +\frac{{G_{\gamma }}_C}{R_e^2}\varphi ,$

(55)$\frac{\partial \theta }{\partial t^{\alpha }}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\frac{1}{P_r}\frac{1}{R_e}\frac{{\partial }^{2}u}{\partial y^2}+\frac{N_b}{R_e}\frac{\partial \theta }{\partial y}\frac{\partial \varphi }{\partial y}+\frac{N_t}{R_e}{\left(\frac{\partial \theta }{\partial y}\right)}^2+\frac{\epsilon }{x}\frac{1}{P_r}\left(A^{*}u+B^{*}\theta \right),$

(56)$\frac{\partial \varphi }{\partial t^{\alpha }}+u\frac{\partial \varphi }{\partial x}+v\frac{\partial \varphi }{\partial y}=\frac{1}{S_c}\frac{1}{R_e}\frac{{\partial }^2\varphi }{\partial y^2}+\frac{N_t}{N_b}\frac{1}{S_c{R}_e}\frac{{\partial }^2\theta }{\partial y^2}+\frac{S_r}{R_e}\frac{{\partial }^2\theta }{\partial y^2}-\gamma \varphi .$

Subject to the same boundary conditions of the classical model.

Justification for Time Fractal Model: The incorporation of the time fractal model (Equations (54)–(56)) in our study is driven by the need to investigate the possible advantages of fractal derivatives in capturing complex temporal patterns in the flow of non-Newtonian fluids across a plate. The parameter $\alpha$ in the time fractional derivatives signifies the fractal time order, which indicates the irregular nature of the temporal progression of fluid dynamics. Fractal derivatives in time provide a more detailed explanation of time-dependent events, encompassing long-term memory effects and non-local dependencies that classical time derivatives may not fully encompass.

Physical meaning of $\alpha$: The significance of $\alpha$ lies in its impact on the temporal progression of the system. A value of α ranging from 0 to 1 suggests subdiffusion, but $\alpha >1$ denotes superdiffusion. The choice of α is contingent upon the particular attributes of the fluid flow and the intended degree of temporal complexity in the model. By incorporating fractal time derivatives, our objective is to investigate if this method may offer a more precise depiction of the dynamics of non-Newtonian fluid flow.

Comparison with Classical Model (42)–(45): To clarify, the classical model (Equations (42)–(45)) describes the non-Newtonian fluid flow at a baseline, and the proposed fractal model is an effort to improve the classical model’s ability to capture complex temporal features by introducing fractal derivatives in time. By adding fractional time derivatives, the fractal model takes the classical model to a new degree of complexity that could be useful in some fluid dynamics situations.

Numerical scheme based on Taylor series: The numerical approaches used in this study are based on Taylor series expansions. However, what sets them apart is their ability to handle fractal temporal derivatives. The numerical solution of fractal differential equations presents distinct difficulties because of fractional derivatives’ non-local and memory-dependent characteristics. We employ Taylor series expansions as a methodological approach to tackle these issues and computationally solve the suggested fractal model.

The novelty of our study is that it investigates fractal time derivatives to improve the temporal precision of simulations of non-Newtonian fluid flow. The suggested model aims to expand the limits of conventional fluid dynamics modelling and explore the advantages of integrating fractal notions into the temporal domain.

In this research, our parametric analysis is useful for several reasons:

• Understanding System Behavior: Our objective is to systematically examine the effects of changing various parameters on the profiles of velocity, temperature, and nanoparticle volume percentage in the boundary layer flow. These parameters include the Prandtl number, Eyring viscosity parameter, and heat generation parameters. Important for comprehending complicated fluid dynamics events, this sheds light on how alterations to these parameters impact the general behavior of the fluid’s flow and thermal properties.

• Validation of Computational Scheme: Parametric analysis validates our computational scheme. We may evaluate our computational approach’s correctness and dependability by comparing our numerical results to theoretical data for specific parameter values. Validating our numerical simulations and findings is crucial to their credibility and robustness.

• Practical Implications and Engineering Applications: Our parametric study lays the framework for understanding non-Newtonian nanofluid behavior under different conditions, but it does not explicitly address practical applications. Thermal management, nanofluid manufacturing, and renewable energy engineers and researchers benefit from such knowledge. Our study helps optimize and develop nanofluid systems by revealing how parameters affect fluid flow and heat transfer.

5. Results and Discussions

A proposed fractal scheme can handle fractal and classical partial differential equations. The scheme is a three-stage scheme, and all stages of the scheme are explicit. It does not require any other scheme to solve difference equations that arise when discretizing the given partial differential equations by the proposed scheme. So, this is one of the advantages of using explicit schemes over implicit schemes. The technique is applicable to time fractal partial differential equations, but not to fractal space partial differential equations. This is because the scheme uses a second-order central spatial discretization, which can be improved by employing a fourth or higher-order spatial discretization. The scheme’s stability analysis is provided, which is applicable to scalar partial differential equations. The stability of the suggested method is assessed using Fourier series analysis. The system exhibits third-order accuracy in fractal time and second-order accuracy in space, thus establishing its consistency. According to the Lax equivalence theorem, it can be verified that it is convergent.

A new fractal computing technique that can solve classical and fractal PDEs is presented in our study. The explicit nature and capability to handle fractal time-dependent boundary layer flow of non-Newtonian nanofluids with space and temperature-dependent heat generation are two areas where this scheme excels over current numerical methods. Its fractal time accuracy is third-order, and its space accuracy is second-order.

The numerical scheme is capable of obtaining the solution of the classical model. But it is used in this work for solving fractal partial differential equations. The effect of dimensionless parameters on velocity, temperature, and nanoparticle volume fraction is similar to those of the classical model. If the value of parameter $\mathsf{\alpha }$ is one, then the classical results can be obtained. As explained earlier, the time fractal derivative can just shift the solution because only the run or denominator of the classical derivative will be changed. So, due to the consideration of fractal derivative, the solution of the classical model will be shifted along the $x$-axis. The effect of fractal derivative can be observed for all involved dimensionless parameters. The physical reason will remain the same as for the solution of the classical model. Still, variation of each parameter, the solution of the classical model will be shifted along the axis. So, physical justification will remain valid for the fractal model.

Figure 1 illustrates the impact of parameter I on the velocity profile. The velocity profile increases as parameter I increases. The velocity profile is enhanced since the diffusion process rises as parameter I grows. Figure 2 displays the effect of the Hartman number on the velocity profile. The velocity profile declines as the Hartman number grows. The decline of the velocity profile is the consequence of the increment of the Lorentz force that resists the velocity of the flow. Figure 3 illustrates the change in the thermal Grashof number across the velocity profile. The velocity profile increases as the thermal Grashof number increases. Since the temperature gradient is one of the factors that drive the flow of mixed convective flows, by enhancing the Grashof number, the temperature gradient rises and increases the flow velocity. Figure 4 illustrates the changes in the thermophoresis parameter across the temperature profile. As the thermophoresis parameter increases, the temperature profile also increases. The rise in the temperature profile increases the cycling of hot and cold particles. As the thermophoresis parameter increases, colder particles move closer to the plate, while hotter and lighter particles move away from the plate towards their surroundings. This movement of particles elevates the temperature of the fluid.

The impact of the Prandtl number on temperature profile is depicted in Figure 5. Temperature profile declines by raising the Prandtl number. The Prandtl number is inversely proportional to thermal diffusivity. A higher Prandtl number causes a lower thermal diffusivity, which in turn leads to a lower thermal conductivity. Consequently, the temperature profile experiences a reduction. The effect of the Soret number on nanoparticle volume fraction is displayed in Figure 6. The nanoparticle volume fraction grows by raising the Soret number. By increasing the Soret number, the nanoparticle volume fraction difference may increase, leading to a rise in nanoparticle volume fraction. Figure 7 shows the Brownian motion parameter on nanoparticle volume fraction. Nanoparticle volume fraction declines by raising the Brownian motion parameter. The difference in nanoparticle volume fraction might increase by incrementing the Brownian motion parameter so that this difference will decrease the nanoparticle volume fraction.

The effect of parameter I and parameter II on the skin friction coefficient is displayed in Figure 8. The skin friction coefficient has both increasing and decreasing behavior by raising parameter I. The skin friction coefficient declines by enhancing parameter II. The effect of the heat source parameter for the space-dependent term and Prandtl number on the local Nusselt number is shown in Figure 9. The Local Nusselt number rises by incrementing the Prandtl number and declines by a raising the coefficient of the dependent term of the heat source. Since the temperature profile rises by raising the heat source parameters and leads to a decrease in heat flux, the local Nusselt number declines. Also, the temperature profile declines by raising the Prandtl number, resulting in heat flux, and so the local Nusselt number rises. Figure 10 shows the influence of the Schmidt number and reaction rate parameter on the local Sherwood number. The local Sherwood number rises as the Schmidt number and reaction rate parameter grows. Since the nanoparticle volume fraction decays by raising the Schmidt number and the reaction rate parameter, this results in growth in mass flux and therefore increases nanoparticle volume fraction.

Figure 11, Figure 12 and Figure 13 show the contour plots for the horizontal velocity component and the flow’s temperature profile over an oscillatory sheet. The effect of the oscillatory boundary can be seen in Figure 11, which shows a contour plot over time and spatial coordinates. Since the oscillatory boundary changes by varying time, its effect is clear along time coordinates. Figure 14 and Figure 15 show the effect of different fractal derivatives on local Nusselt numbers and local Sherwood numbers with different values of $\alpha ,$ respectively.

Table 2 compares three numerical schemes for solving the two-dimensional parabolic equation, which appeared as the first example in [49]. The suggested scheme exhibits third-order accuracy in the temporal dimension. The forward Euler scheme is the first-order scheme, while the second-order Runge–Kutta scheme is the third. The spatial discretization is carried out using the second-order central difference formula. This comparison demonstrates that the proposed approach yields a lower margin of error compared to the other two alternatives.

Building on our previous discussion of the data, we would want to draw attention to a number of important discoveries and what they mean.

Parametric Analysis and Physical Insights: A comprehensive parametric analysis examined how dimensionless parameters affected velocity, temperature, and nanoparticle volume fraction profiles. Our findings show the complex relationship between these parameters and fluid flow and thermal behavior. The effect of parameters like the Hartman number on velocity profiles shows how magnetic forces resist fluid motion, while thermal Grashof and Prandtl numbers show how temperature gradients and thermal diffusivity affect heat transfer processes.

Comparative Analysis and Validation: The forward Euler and second-order Runge–Kutta systems were compared to our fractal scheme. Our results show that the suggested method improves accuracy and reduces errors for two-dimensional parabolic equations. Our computational technique is validated and emphasized as superior in terms of accuracy and efficiency in this comparative analysis, which advances the state-of-the-art in numerical simulations of complicated fluid dynamics issues.

Integration with State-of-the-Art Research: We will broaden our discussion to incorporate references to current studies that examine similar phenomena or use comparable numerical approaches, allowing for a better contextualization of our findings within the larger framework of fluid dynamics research. Our findings can be better understood in relation to the current state of the art, and our contributions can be more prominently displayed if we compare them to previous literature.

6. Conclusions

This paper studies an explicit computational approach to solve the mixed convective flow of Prandtl-Eyring nanofluid with the effect of space and temperature-dependent heat generation. The proposed three-stage scheme, constructed based on the fractal Taylor series, has proven effective in solving the fractal time-dependent partial differential equations. It is vital to determine the developed scheme’s stability and usefulness. A computational scheme for non-Newtonian Prandtl-Eyring nanofluid flow over solid surfaces has been proposed. Its convergence is also provided for the fractal time-dependent partial differential equations system. The proposed scheme has been employed for fractal time terms, and second-order central difference discretization is used for space discretization. The graphs illustrate the impact of various factors on velocity, temperature, and nanoparticle volume fraction. The concluding points can be stated in the following way:

• The velocity profile was enhanced as the parameter I grew.

• The rising coefficient of the space-dependent heat generation term increased the temperature profile.

• The rising Brownian motion parameter declined the nanoparticle volume fraction.

• Local Sherwood number increased with Schmidt number and response rate parameter.

By providing a new way to deal with the problems caused by mixed convective flows of nanofluids with non-Newtonian properties, this study adds to the larger area of computational fluid dynamics. The results are more relevant to real-world situations since space and temperature-dependent heat generation are included in the simulations, which makes them more realistic.

Future Directions: A potential future research area could involve expanding the proposed scheme to explore additional sophisticated fluid dynamics scenarios while considering additional parameters and real-world boundary conditions. Extending the proposed scheme to investigate other complex fluid dynamics scenarios, taking into account extra parameters and real-world boundary conditions, could be a future study direction. This study only scratches the surface of the complex parameter interaction within the system; future research may shed light on the system’s behavior and point the way toward better predictive models.

This paper offers a computational framework [50,51,52] to tackle a particular issue in nanofluid mechanics and contributes to the broader field of numerical methods in fluid dynamics. The discoveries provided here provide a basis for future efforts to understand fluid flow’s details in many practical settings, which are crucial for progress in engineering and technology.

Author Contributions

Conceptualization, methodology, and analysis, Y.N.; funding acquisition, A.S.B.; investigation, M.M.; methodology, Y.N.; project administration, K.A.; resources, K.A.; supervision, M.S.A.; visualization, M.M.; writing—review and editing, M.S.A.; proofreading and editing, A.S.B. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The manuscript included all required data and implementing information.

Acknowledgments

The authors wish to express their gratitude to Prince Sultan University for facilitating the publication of this article through the Theoretical and Applied Sciences Lab.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Effect of parameter I on velocity profile using [Forumla omitted. See PDF.].

Figure 2. Effect of Hartmann number on velocity profile using [Forumla omitted. See PDF.].

Figure 3. Effect of thermal Grashof number on velocity profile using [Forumla omitted. See PDF.].

Figure 4. Effect of thermophoresis parameter on temperature profile using [Forumla omitted. See PDF.].

Figure 5. Effect of Prandtl number on temperature profile using [Forumla omitted. See PDF.].

Figure 6. Effect of Soret number on nanoparticle volume fraction profile using [Forumla omitted. See PDF.].

Figure 7. Effect of Brownian motion parameter on nanoparticle volume fraction profile using [Forumla omitted. See PDF.].

Figure 8. Effect of parameter I and parameter II on skin friction coefficient using [Forumla omitted. See PDF.].

Figure 9. Effect of coefficient of space-dependent term of heat generation and Prandtl number on local Nusselt number using [Forumla omitted. See PDF.].

Figure 10. Effect of Schmidt number and reaction rate parameter on local Sherwood number using [Forumla omitted. See PDF.].

Figure 11. Contour plot for the horizontal component of velocity along time and spatial coordinates using [Forumla omitted. See PDF.].

Figure 12. Contour plot for the horizontal component of velocity along spatial coordinates using [Forumla omitted. See PDF.].

Figure 13. Contour plot for temperature profile along spatial coordinates using [Forumla omitted. See PDF.].

Figure 14. Variation of fractal derivative parameter on local Nusselt number using [Forumla omitted. See PDF.].

Figure 15. Variation of fractal derivative parameter on local Sherwood number using [Forumla omitted. See PDF.].

Table 1

Fractional derivatives have been studied from numerous aspects in boundary layer problems.

Fluid	Type of Derivative	Type of Magnetic Field	Type of Solution	Researchers	Case Study
Single nanofluid (Water-CuWater-Ag Water-$A{l}_2{O}_3$ Water-$Ti{O}_2$)	Spatial fractional	Without magnetic field	Numerical (finite difference)	Pan et al. [35]	Boundary layer flow in a porous media
Second-grade fluid	Time fractional	Inclined magnetic field	Numerical (Laplace along with Zakian’s algorithm)	Tassaddiq [36]	Boundary layer flow along an inclined heated plate
Viscoelastic fluid	Time fractional	Vertical magnetic field	Numerical (finite difference)	Chen et al. [37]	Boundary layer flow over a stretching sheet
Maxwell fluid	Time fractional	Without magnetic field	Numerical (finite difference)	Yang et al. [38]	Stretching sheet with variable thickness
Viscoelastic fluid	Spatial fractional	Without magnetic field	Numerical (finite difference)	Li et al. [39]	Boundary layer over a permeable surface
Sisko nanofluid	Time fractional	Without magnetic field	Numerical (finite difference)	Shen et al. [40]	Boundary layer flow over a continuously moving plate
Maxwell fluid	Time fractional	Without magnetic field	Numerical (finite difference)	Liu et al. [41]	Boundary layer over a moving plate
Single nanofluid (Water-SWCNTs Water-MWCNTs)	Time fractional	Vertical magnetic field	Numerical (Joint of finite-difference discretization and L1 algorithm)	Anwar et al. [42]	Boundary layer flow induced due to a stretching sheet

Table 2

Comparison of three schemes using ${\mathrm{N}}_{\mathrm{x}}=50={\mathrm{N}}_{\mathrm{y}}(\mathrm{N}\mathrm{o}.\mathrm{o}\mathrm{f} \mathrm{g}\mathrm{r}\mathrm{i}\mathrm{d} \mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s})$ ${\mathrm{t}}_{\mathrm{f}}(\mathrm{F}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l} \mathrm{T}\mathrm{i}\mathrm{m}\mathrm{e})=0.1$.

$∆\mathit{t}$	Proposed		1st Order Scheme		2nd Order Scheme
	${\mathbf{L}}_2$ Error	Conv. Order	${\mathbf{L}}_2$ Error	Conv. Order	${\mathbf{L}}_2$ Error	Conv. Order
$\raisebox{1ex}{$0.1$}\!\left/ \!\raisebox{-1ex}{$1500$}\right.$	4.27 × 10−4	-	6.29 × 10−4	-	5.31 × 10−4	-
$\raisebox{1ex}{$0.1$}\!\left/ \!\raisebox{-1ex}{$1750$}\right.$	4.19 × 10−4	0.1227	5.93 × 10−4	0.3823	5.09 × 10−4	0.2745
$\raisebox{1ex}{$0.1$}\!\left/ \!\raisebox{-1ex}{$2000$}\right.$	4.13 × 10−4	0.1080	5.66 × 10−4	0.3490	4.92 × 10−4	1.2805
$\raisebox{1ex}{$0.1$}\!\left/ \!\raisebox{-1ex}{$2250$}\right.$	4.09 × 10−4	0.0826	5.45 × 10−4	0.3210	4.79 × 10−4	0.2274

References

1. Aziz, T.; Fatima, A.; Khalique, C.M.; Mahomed, F.M. Prandtl’s boundary layer equation for two-dimensional flow: Exact solutions via the simplest equation method. Math. Probl. Eng.; 2013; 2013, 724385. [DOI: https://dx.doi.org/10.1155/2013/724385]

2. Sankad, G.; Ishwar, M.; Dhange, M. Varying wall temperature and thermal radiation effects on MHD boundary layer liquid flow containing gyrotactic microorganisms. Partial. Differ. Equ. Appl. Math.; 2021; 4, 100092.

3. Hussain, M.; Ghaffar, A.; Ali, A.; Shahzad, A.; Nisar, K.S.; Alharthi, M.R.; Jamshed, W. MHD thermal boundary layer flow of a Casson fluid over a penetrable stretching wedge in the existence of nonlinear radiation and convective boundary condition. Alex. Eng. J.; 2021; 60, pp. 5473-5483. [DOI: https://dx.doi.org/10.1016/j.aej.2021.03.042]

4. Abedi, H.; Sarkar, S.; Johansson, H. Numerical modelling of neutral atmospheric boundary layer flow through heterogeneous forest canopies in complex terrain (a case study of a Swedish wind farm). Renew. Energy; 2021; 180, pp. 806-828. [DOI: https://dx.doi.org/10.1016/j.renene.2021.08.036]

5. Yang, S.; Liu, L.; Long, Z.; Feng, L. Unsteady natural convection boundary layer flow and heat transfer past a vertical flat plate with novel constitution models. Appl. Math. Lett.; 2021; 120, 107335. [DOI: https://dx.doi.org/10.1016/j.aml.2021.107335]

6. Long, Z.; Liu, L.; Yang, S.; Feng, L.; Zheng, L. Analysis of Marangoni boundary layer flow and heat transfer with novel constitution relationships. Int. Commun. Heat Mass Transf.; 2021; 127, 105523. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2021.105523]

7. Zhang, X.; Tang, Y.; Zhang, F.; Lee, C.S. A novel aluminum—Graphite dual-ion battery. Adv. Energy Mater.; 2016; 6, 1502588. [DOI: https://dx.doi.org/10.1002/aenm.201502588]

8. Hosseini, S.M.; Safaei, M.R.; Goodarzi, M.; Alrashed, A.A.; Nguyen, T.K. New temperature, interfacial shell dependent dimensionless model for thermal conductivity of nanofluids. Int. J. Heat Mass Transf.; 2017; 114, pp. 207-210. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.06.061]

9. Ahmadi, M.H.; Mohseni-Gharyehsafa, B.; Ghazvini, M.; Goodarzi, M.; Jilte, R.D.; Kumar, R. Comparing various machine learning approaches in modeling the dynamic viscosity of CuO/water nanofluid. J. Therm. Anal. Calorim.; 2020; 139, pp. 2585-2599. [DOI: https://dx.doi.org/10.1007/s10973-019-08762-z]

10. Bahiraei, M.; Salmi, H.K.; Safaei, M.R. Effect of employing a new biological nanofluid containing functionalized graphene nanoplatelets on thermal and hydraulic characteristics of a spiral heat exchanger. Energy Convers. Manag.; 2019; 180, pp. 72-82. [DOI: https://dx.doi.org/10.1016/j.enconman.2018.10.098]

11. Waqas, H.; Farooq, U.; Khan, S.A.; Alshehri, H.M.; Goodarzi, M. Numerical analysis of dual variable of conductivity in bioconvection flow of Carreau–Yasuda nanofluid containing gyrotactic motile microorganisms over a porous medium. J. Therm. Anal. Calorim.; 2021; 145, pp. 2033-2044. [DOI: https://dx.doi.org/10.1007/s10973-021-10859-3]

12. Wang, X.; Li, C.; Zhang, Y.; Said, Z.; Debnath, S.; Sharma, S.; Yang, M.; Gao, T. Influence of texture shape and arrangement on nanofluid minimum quantity lubrication turning. Int. J. Adv. Manuf. Technol.; 2021; 119, pp. 631-646. [DOI: https://dx.doi.org/10.1007/s00170-021-08235-4]

13. Xie, Y.; Meng, X.; Mao, D.; Qin, Z.; Wan, L.; Huang, Y. Homogeneously dispersed graphene nanoplatelets as long-term corrosion inhibitors for aluminum matrix composites. ACS Appl. Mater. Interfaces; 2021; 13, pp. 32161-32174. [DOI: https://dx.doi.org/10.1021/acsami.1c07148] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34185498]

14. Abu-Hamdeh, N.H.; Alsulami, R.A.; Rawa, M.J.; Alazwari, M.A.; Goodarzi, M.; Safaei, M.R. A Significant Solar Energy Note on Powell-Eyring Nanofluid with Thermal Jump Conditions: Implementing Cattaneo-Christov Heat Flux Model. Mathematics; 2021; 9, 2669. [DOI: https://dx.doi.org/10.3390/math9212669]

15. Sajid, T.; Jamshed, W.; Shahzad, F.; Eid, M.R.; Alshehri, H.M.; Goodarzi, M.; Akgül, E.K.; Nisar, K.S. Micropolar fluid past a convectively heated surface embedded with nth order chemical reaction and heat source/sink. Phys. Scr.; 2021; 96, 104010. [DOI: https://dx.doi.org/10.1088/1402-4896/ac0f3e]

16. Maleki, H.; Safaei, M.R.; Togun, H.; Dahari, M. Heat transfer and fluid flow of pseudo-plastic nanofluid over a moving permeable plate with viscous dissipation and heat absorption/generation. J. Therm. Anal. Calorim.; 2019; 135, pp. 1643-1654. [DOI: https://dx.doi.org/10.1007/s10973-018-7559-2]

17. Shankar, U.; Naduvinamani, N. Magnetized squeezed flow of time-dependent Prandtl-Eyring fluid past a sensor surface. Heat Transf. Asian Res.; 2019; 48, pp. 2237-2261. [DOI: https://dx.doi.org/10.1002/htj.21482]

18. Al-Kaabi, W.; Al-Khafajy, D.G.S. Radiation and Mass Transfer Effects on Inclined MHD Oscillatory Flow for Prandtl-Eyring Fluid through a Porous Channel. Al-Qadisiyah J. Pure Sci.; 2021; 26, pp. 347-363.

19. Hayat, T.; Ullah, I.; Muhammad, K.; Alsaedi, A. Gyrotactic microorganism and bio-convection during flow of Prandtl-Eyring nanomaterial. Nonlinear Eng.; 2021; 10, pp. 201-212. [DOI: https://dx.doi.org/10.1515/nleng-2021-0015]

20. Waqas, H.; Farooq, U.; Alshehri, H.M.; Goodarzi, M. Marangoni-bioconvectional flow of Reiner–Philippoff nanofluid with melting phenomenon and nonuniform heat source/sink in the presence of a swimming microorganisms. Math. Methods Appl. Sci.; 2021; [DOI: https://dx.doi.org/10.1002/mma.7727]

21. Alazwari, M.A.; Abu-Hamdeh, N.H.; Goodarzi, M. Entropy Optimization of First-Grade Viscoelastic Nanofluid Flow over a Stretching Sheet by Using Classical Keller-Box Scheme. Mathematics; 2021; 9, 2563. [DOI: https://dx.doi.org/10.3390/math9202563]

22. Munjam, S.R.; Gangadhar, K.; Seshadri, R.; Rajeswar, M. Novel technique MDDIM solutions of MHD flow and radiative Prandtl-Eyring fluid over a stretching sheet with convective heating. Int. J. Ambient. Energy; 2021; 43, pp. 4850-4859. [DOI: https://dx.doi.org/10.1080/01430750.2021.1922498]

23. Jamshed, W.; Kumar, V.; Kumar, V. Computational examination of Casson nanofluid due to a nonlinear stretching sheet subjected to particle shape factor: Tiwari and Das model. Numer. Methods Partial. Differ. Equ.; 2020; 38, pp. 848-875. [DOI: https://dx.doi.org/10.1002/num.22705]

24. Jamshed, W.; Mishra, S.; Pattnaik, P.; Nisar, K.S.; Devi, S.S.U.; Prakash, M.; Shahzad, F.; Hussain, M.; Vijayakumar, V. Features of entropy optimization on viscous second grade nanofluid streamed with thermal radiation: A Tiwari and Das model. Case Stud. Therm. Eng.; 2021; 27, 101291. [DOI: https://dx.doi.org/10.1016/j.csite.2021.101291]

25. Abdelmalek, Z.; Hussain, A.; Bilal, S.; Sherif, E.-S.M.; Thounthong, P. Brownian motion and thermophoretic diffusion influence on thermophysical aspects of electrically conducting viscoinelastic nanofluid flow over a stretched surface. J. Mater. Res. Technol.; 2020; 9, pp. 11948-11957. [DOI: https://dx.doi.org/10.1016/j.jmrt.2020.08.014]

26. Pan, M.; Zheng, L.; Liu, F.; Zhang, X. Modeling heat transport in nanofluids with stagnation point flow using fractional calculus. Appl. Math. Model.; 2016; 40, pp. 8974-8984. [DOI: https://dx.doi.org/10.1016/j.apm.2016.05.044]

27. Li, S.N.; Cao, B.Y. Fractional Boltzmann transport equation for anomalous heat transport and divergent thermal conductivity. Int. J. Heat Mass Transf.; 2019; 137, pp. 84-89. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.03.120]

28. Santina, D.; Asif, M.; Aljawi, S.; Mlaiki, N. Application of the inverse Laplace transform techniques to solve the generalized Bagley–Torvik equation including Caputo’s fractional derivative. Partial. Differ. Equ. Appl. Math.; 2024; 10, 100664. [DOI: https://dx.doi.org/10.1016/j.padiff.2024.100664]

29. Chen, X.; Ye, Y.; Zhang, X.; Zheng, L. Lie-group similarity solution and analysis for fractional viscoelastic MHD fluid over a stretching sheet. Comput. Math. Appl.; 2018; 75, pp. 3002-3011. [DOI: https://dx.doi.org/10.1016/j.camwa.2018.01.028]

30. Caputo, M. Linear models of dissipation whose Q is almost frequency independent. Ann. Geophys.; 1966; 19, pp. 383-393. [DOI: https://dx.doi.org/10.4401/ag-5051]

31. Chakrabarty, A.K.; Roshid, M.M.; Rahaman, M.M.; Abdeljawad, T.; Osman, M.S. Dynamical analysis of optical soliton solutions for CGL equation with Kerr law nonlinearity in classical, truncated M-fractional derivative, beta fractional derivative, and conformable fractional derivative types. Results Phys.; 2024; 60, 107636. [DOI: https://dx.doi.org/10.1016/j.rinp.2024.107636]

32. Asjad, M.I.; Sarwar, N.; Hafeez, M.B.; Sumelka, W.; Muhammad, T. Advancement of non-newtonian fluid with hybrid nanoparticles in a convective channel and prabhakar’s fractional derivative—Analytical solution. Fractal Fract.; 2021; 5, 99. [DOI: https://dx.doi.org/10.3390/fractalfract5030099]

33. Asjad, M.I.; Ali, R.; Iqbal, A.; Muhammad, T.; Chu, Y.M. Application of water based drilling clay-nanoparticles in heat transfer of fractional Maxwell fluid over an infinite flat surface. Sci. Rep.; 2021; 11, 25. [DOI: https://dx.doi.org/10.1038/s41598-021-98066-w]

34. Asjad, M.I.; Basit, A.; Iqbal, A.; Shah, N.A. Advances in transport phenomena with nanoparticles and generalized thermal process for vertical plate. Phys. Scr.; 2021; 96, 15. [DOI: https://dx.doi.org/10.1088/1402-4896/ac102b]

35. Pan, M.; Zheng, L.; Liu, F.; Liu, C.; Chen, X. A spatial-fractional thermal transport model for nanofluid in porous media. Appl. Math. Model.; 2018; 53, pp. 622-634. [DOI: https://dx.doi.org/10.1016/j.apm.2017.08.026]

36. Tassaddiq, A. MHD flow of a fractional second grade fluid over an inclined heated plate. Chaos Solitons Fractals; 2019; 123, pp. 341-346. [DOI: https://dx.doi.org/10.1016/j.chaos.2019.04.029]

37. Chen, X.; Yang, W.; Zhang, X.; Liu, F. Unsteady boundary layer flow of viscoelastic MHD fluid with a double fractional Maxwell model. Appl. Math. Lett.; 2019; 95, pp. 143-149. [DOI: https://dx.doi.org/10.1016/j.aml.2019.03.036]

38. Yang, W.; Chen, X.; Zhang, X.; Zheng, L.; Liu, F. Flow and heat transfer of double fractional Maxwell fluids over a stretching sheet with variable thickness. Appl. Math. Model.; 2020; 80, pp. 204-216. [DOI: https://dx.doi.org/10.1016/j.apm.2019.11.017]

39. Li, B.; Liu, F. Boundary layer flows of viscoelastic fluids over a non-uniform permeable surface. Comput. Math. Appl.; 2020; 79, pp. 2376-2387. [DOI: https://dx.doi.org/10.1016/j.camwa.2019.11.003]

40. Shen, M.; Chen, L.; Zhang, M.; Liu, F. A renovated Buongiorno’s model for unsteady Sisko nanofluid with fractional Cattaneo heat flux. Int. J. Heat Mass Transf.; 2018; 126, pp. 277-286. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.05.131]

41. Liu, L.; Feng, L.; Xu, Q.; Zheng, L.; Liu, F. Flow and heat transfer of generalized Maxwell fluid over a moving plate with distributed order time fractional constitutive models. Int. Commun. Heat Mass Transf.; 2020; 116, 415. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2020.104679]

42. Anwar, M.S.; Ahmad, R.T.M.; Shahzad, T.; Irfan, M.; Ashraf, M.Z. Electrified fractional nanofluid flow with suspended carbon nanotubes. Comput. Math. Appl.; 2020; 80, pp. 1375-1386. [DOI: https://dx.doi.org/10.1016/j.camwa.2020.07.005]

43. Sun, Z.Z.; Wu, X. A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math.; 2006; 56, pp. 193-209. [DOI: https://dx.doi.org/10.1016/j.apnum.2005.03.003]

44. She, M.; Li, D.; Sun, H.W. A transformed L1 method for solving the multi-term time-fractional diffusion problem. Math. Comput. Simul.; 2022; 193, pp. 584-606. [DOI: https://dx.doi.org/10.1016/j.matcom.2021.11.005]

45. Yuan, W.; Li, D.; Zhang, C. Linearized Transformed L1 Galerkin FEMs with Unconditional Convergence for Nonlinear Time Fractional Schrödinger Equations. Numer. Math. Theory Methods Appl.; 2023; 16, pp. 348-369. [DOI: https://dx.doi.org/10.4208/nmtma.OA-2022-0087]

46. Lu, Z.; Fan, W. A fast algorithm for multi-term time-space fractional diffusion equation with fractional boundary condition. Numer. Algorithms; 2024; pp. 1-24. [DOI: https://dx.doi.org/10.1007/s11075-024-01830-y]

47. Chen, W.; Sun, H.; Zhang, X.; Korošak, D. Anomalous diffusion modeling by fractal and fractional derivatives. Comput. Math. Appl.; 2010; 59, pp. 1754-1758. [DOI: https://dx.doi.org/10.1016/j.camwa.2009.08.020]

48. Khan, Z.A.; Shah, K.; Abdalla, B.; Abdeljawad, T. A numerical Study of Complex Dynamics of a Chemostat Model Under Fractal-Fractional Derivative. Fractals; 2023; 31, 2340181. [DOI: https://dx.doi.org/10.1142/S0218348X23401813]

49. Li, F.-l.; Wu, Z.-k.; Ye, C.-r. A finite difference solution to a two-dimensional parabolic inverse problem. Appl. Math. Model.; 2012; 36, pp. 2303-2313. [DOI: https://dx.doi.org/10.1016/j.apm.2011.08.025]

50. Arif, M.S.; Abodayeh, K.; Nawaz, Y. A finite difference explicit-implicit scheme for fractal heat and mass transportation of Williamson nanofluid flow in quantum calculus. Numer. Heat Transf. Part A Appl.; 2024; pp. 1-23. [DOI: https://dx.doi.org/10.1080/10407782.2024.2308753]

51. Arif, M.S.; Abodayeh, K.; Nawaz, Y. The modified finite element method for heat and mass transfer of unsteady reacting flow with mixed convection. Front. Phys.; 2022; 10, 952787. [DOI: https://dx.doi.org/10.3389/fphy.2022.952787]

52. Nawaz, Y.; Arif, M.S.; Abodayeh, K. An explicit-implicit numerical scheme for time fractional boundary layer flows. Int. J. Numer. Methods Fluids; 2022; 94, pp. 920-940. [DOI: https://dx.doi.org/10.1002/fld.5078]