# Finite difference schemes for MHD mixed convective Darcy–forchheimer flow of Non-Newtonian fluid over oscillatory sheet: A computational study

^{1}Department of Mathematics, Air University, Islamabad, Pakistan^{2}Department of Mathematics and Sciences, College of Humanities and Sciences, Prince Sultan University, Riyadh, Saudi Arabia^{3}Department of Mathematics, University of Wah, Wah Cant, Pakistan

This contribution proposes two third-order numerical schemes for solving time-dependent linear and non-linear partial differential equations (PDEs). For spatial discretization, a compact fourth-order scheme is deliberated. The stability of the proposed scheme is set for scalar partial differential equation, whereas its convergence is specified for a system of parabolic equations. The scheme is applied to linear scalar partial differential equation and non-linear systems of time-dependent partial differential equations. The non-linear system comprises a set of governing equations for the heat and mass transfer of magnetohydrodynamics (MHD) mixed convective Casson nanofluid flow across the oscillatory sheet with the Darcy–Forchheimer model, joule heating, viscous dissipation, and chemical reaction. It is noted that the concentration profile is escalated by mounting the thermophoresis parameter. Also, the proposed scheme converges faster than the existing Crank-Nicolson scheme. The findings that were provided in this study have the potential to serve as a helpful guide for investigations into fluid flow in closed-off industrial settings in the future.

## 1 Introduction

Electromagnetism and fluid mechanics share a modern scientific field called magnetohydrodynamics (MHD). The electromagnetic force, which has been widely exploited in many industries, such as semiconductor crystal growth, liquid metal blankets in nuclear fusion reactors, and steel-making processes, significantly hinders the flow of conducting materials.

Levitation of drops, dynamo modelling of planets, electromagnetic pumps, and other technologies. In light of these advancements, it is crucial to research significant MHD flows such as rotating, free-surface, and natural convection in rectangular ducts. Numerical analyses are being utilized to examine MHD flows that are more composite due to recent advancements in computational approaches. This Special Issue emphasized numerical and computational methods for deciphering and analyzing more complex MHD processes [1]. Miura [2] examined magneto hydrodynamic incompressible turbulence driven by gyro-viscous terms and Hall by numerical simulations of freely decaying, nearly isotropic, and homogeneous turbulent flow. He discovered that the Hall and gyro-viscous variables modify the energy transfer in the equations of motion to be predominantly forward transfer, while the magnetic energy transfer stays dominating.

The problem to be tackled in the present work is the Casson nanofluid flow across oscillatory surfaces with effects of joule heating, thermal radiations, and chemical reaction. Casson nan fluid is preferable to Newtonian-based nan fluid flow as cooling and friction-reducing agents [3]. Casson introduced the Casson fluid model for the movement of viscoelastic liquids in 1959. Stress at which no flow occurs should be provided by the Casson fluid, a shear-thinning fluid, which should have infinite viscosity at zero rates of shear and zero viscosity at infinite rates of shear. Casson fluid includes honey, jelly, sauce, soup, and other things [4]. Ali et al.'s [5] investigation of the Casson fluid flow on a tilted sheet included the Soret-Dufour effects. Manideep et al. [6] studied the fluid Casson flow on vertically inclined sheets. Shamshuddin et al. [7] conducted a statistical analysis of a chemical reaction’s impact on the Casson fluid flow on a sloped plate [8]; Vijayaragavan and Kavitha estimated the Casson fluid flow over a tilted plate. Taking into account the current hall, Prasad et al. [9] studied Casson fluid flow on an inclined sheet. The inclined Casson fluid flow on a permeable sheet was studied by Jain and Parmar [10]. The angle effect was considered in Sailaja et al. [11].'s study of the Casson fluid flow on a vertical sheet. Casson fluid flow was investigated by Rawi et al. [12] on a slanted sheet to take nanoparticles into account. Rauju et al. [13] discussed Casson fluid flow on a vertically inclined sheet in their study. The Casson fluid model performs better with blood.

According to the literature, chemical reactions can be divided into boundary or single-phase volumetric reactions, while homo-heterogeneous reactions depend on both. In cases like the creation of smog, first-order homogeneous chemical processes are discussed. The hydrometallurgical and chemical industries are particularly affected by the results of chemical processes involving mass transfer. Since they are inversely correlated with species concentration, reaction rates are frequently thought of as first order. Outside chemicals trigger chemical reactions in the water or the atmosphere. Heat is frequently produced between two species when chemical processes occur, the existence of reactions, and viscous dissipation. The unsteady MHD flow over an extended surface with varying viscosity was studied by Kishan and Hunegnaw [14] while it was embedded in an incredibly porous material. By involving energy and reaction over a parabolic surface using a finite element approach, Y. M. Chu et al. [15] examine the enhancement of thermal energy and solute particles using hybrid nanoparticles. References [16, 17] provide additional details.

A phenomenon of warmth transmission by radiation could be thermal radiation. Due to the high thermal effects involved in the production of electrical engineering and industrial significance can be seen in the cooling of metallic items, paper plates, petroleum pumps, and chips and their effects on MHD [3, 18–20] flows. Applications of radiation effects in physics and engineering processes are considerable. It significantly affects space technology, which studies the thermal effects of various fluxes, and processes that warm solids or liquids. The impacts of thermal radiation might be extremely important in controlling heat transfer in the polymer manufacturing industry, where the value of the materials is focused on heat-maintaining qualities. A bi-viscosity fluid confined in a highly trapezoidal cavity is analyzed by Nasir et al. [21] under the influence of thermal and magnetic forces. According to Kumar et al., the non-linear thermal radiation influence needs to be studied to analyze MHD shear-thickening fluid motion through stretchy geometry, according to Kumar et al. [22]. The effects of thermal radiations on heat transfer for Casson fluid flow in porous media across a non-linearly growing surface are examined by Pramanik [23]. Refs [24–27] provide more details about recent investigations.

Joule heating is when current energy is modified to heat because of flow in resistance. It's a good range of applications, like in industries and technology. The Joule heating effect occurs when an electric current flows through a conductor and causes heat to be generated. The Joule heating effect is used in many industrial processes and devices. In addition to food processing and scanning microscopes, infrared-thermal images, electrical radiative heaters, light bulbs, laboratory water baths, electric tabletop hotplates, clothes irons, hand tools, fan heaters, hair dryers, and cartridge heaters are among the several uses for heating effects of a joule. The researchers and engineers reviewed the Joule heating effect in their respective fields of study in light of the abovementioned applications. The fluid temperature rose for the joule heating parameter for a hyperbolic tangent nanofluid in porous media in a study by Alaidrous et al. [28]. In addition to taking into account the features of the Cattaneo-Christove heat flux theory, Hafeez et al. [29] explored the Oldroyd-B nan non-Newtonian fluid flow. According to this work, the concentration of Oldroyd-B nanofluid increases for thermophoresis. An investigation by Khan et al. [30] determined that the force field parameter increases the nan fluid velocity profile in Casson nanofluid flow through a rotating cylinder that is hydromagnetic in three dimensions, with consequences of entropy and Joule heating. Kumar et al. [31] used the shot method and the RK4 scheme to examine the effects of the Joule heating effect on the flow of Williamson nano fluid toward the stretching sheet. The entropy of the system rises as the Brinkman parameter is advanced. As a cylinder shrinks, a highly hybrid nanofluid flows through it. Khashi’ie et al. [32] investigated the effects of Joule heating and warmth transmission. Uddin et al. [33] identified the multiple features of the Joule heating effect on the mixed convection MHD flow. Our investigation leads us to the conclusion that the fabric magnetic parameter enhanced the skin friction coefficient. Kazemi et al. [34] proposed an analytical solution for the fourth-grade nan fluid flow along the duct walls in viscous dissipation and Joule heating forces. According to their research, heat transport *via* duct walls may reduce by up to 40% when the Hartman number is raised. Rasheed et al. [35] investigated the Jeffery nano liquid MHD boundary flow over an extended cylinder with heat generation and absorption consequences. The Darcy–Forchheimer flow of Williamson fluid over a Riga plate under the effects of suction/injection has been explained in [36]. The heat and mass transfer were also considered with the effect of heat consumption/generation. For solving the ODEs in the presented phenomenon, both numerical and analytical studies were carried out. The numerical solution was based on applying Matlab solver bvp4c, and the analytical study was based on the analytical method named homotopy analysis method. From the obtained results, it was seen that fluid temperature was raised by increasing radiation parameters. A Buongiorno nanofluid model has been considered in [37] for investigating electro-magneto hydrodynamic flow over a Riga plate. A nanofluid model was considered with Cattaneo-Christov’s heat flux and generalized Fick’s law. A study of MHD micropolar ferrofluid flow over a stretching surface has been given in [38]. Using the similarity transformations, the governing equations of flow were reduced into a set of ODEs and solved by Matlab solver bvp4c. It was observed that stream function and velocity declined by enhancing magnetic parameters. Matrix nanocomposites are materials with various applications because of their thermophysical properties.

An all-encompassing numerical study of the many physical features of steady MHD Von Kármán’s flows of chemically reacting nanofluids that can occur over a rotating disk in the presence of a radially applied magnetic field has been performed in [39]. This is done under the assumption that the disk surface is impermeable and heated convectively, in which case the vertical nanoparticles’ mass flux has practically vanished. Examining how the generalized heat transport affects the free convection flow of a viscous nanofluid in a cylinder is what [40] is all about. The heat transfer equation is generalized using the generalized Atangana–Baleanu time-fractional differential operator.

The study for the dispersion of matric nanocomposite material on magnetized nanofluid flow over the coaxial disks was given in [41]. The prominence of the permeability function was also examined numerically. The calculated results showed that the increasing nanolayer thickness had escalated thermal phenomena and effective nanolayer thermal conductivity. It was also noticed that the Nusselt number was enhanced by incrementing hybrid nanoparticles volume fraction. Ferrofluid is made up of tiny magnetic particles and uses in inertial and viscous damping, dynamic sealing and magnetic drug targeting. The study of this fluid has been given in [38], and the stretching sheet caused the flow. For solving ODEs obtained by applying transformations to the governing equations of the fluid flow, Matlab solver bvp4c has been employed. It was observed that stream function and velocity declined by enhancing magnetic parameters. In [42], an infinite medium was used to investigate a thermoelastic diffusion model of fractional order. An inclined temperature field was transmitted through the spherical cavity. The eigenvalue technique solved the differential equation in the form of vector-matric equation in the Laplace transform domain. An analytical solution for the displacement, temperature, and concentration was found in the Laplace domain. A finite element method is one of the numerical methods that can be considered to solve ordinary and partial differential equations. This method has been considered in [43] to study the dual-phase-lag model on thermoelastic interaction subjected to a ramp type of heating. The method was proposed to solve the problem and found the solutions for displacement, temperature and stress. A comparison for some results was also given with three theories. Some more recent work on MHD flow can be found [44–46].

The motivation of this contribution is based on the modifications of analytical methods which have been done earlier, and this effort is not to modify some numerical method. The modification of the numerical method is the main objective of this work and its application to some applied problems in science or engineering. The numerical technique is applied to Casson fluid flow over the oscillatory plate. The non-Newtonian MHD fluid has applications in heat storage beds, textile industries, irrigation problems and polymer composite industries.

In literature, numerous numerical methods exist for solving different types of problems in science and engineering. Some methods discretize only the time variable in the time-dependent PDEs. The time discretization methods can be classified into two types. One of these types of methods is an explicit class of methods, and the second is an implicit class of methods. Explicit methods mostly have smaller stability regions than those provided by some implicit methods. Explicit methods are easy to employ for any non-linear problem, but implicit types of methods require linearized differential equations. Also, the systems of equations obtained by employing implicit methods are required to solve by any software solver solving equations, or these systems require any other iterative method. In this contribution, the Gauss-Seidel iterative method is employed for solving equations after applying the proposed numerical scheme. The proposed numerical schemes are constructed on three-time levels. The first stage is constructed on two-time levels, whereas the second stage uses the information of three-time levels to find the solution at some particular time level. So, a second stage requires any other method based on two levels to find the solution at the first time level. The scheme is applied to two problems: one is a linear problem, and the second is a system of non-linear PDEs. The system of non-linear parabolic equations is obtained from the phenomenon of nanofluid flow over the oscillatory sheet. Effects of heat and mass transfer, thermal radiation, viscous dissipation, and joule heating are also considered. The system is first reduced into dimensionless PDEs and later solved by the proposed scheme with fourth-order compact discretization in space. Analytical and numerical methods can be used to solve considered flow problems. The numerical methods may comprise the finite element method and high-order finite difference approaches that can be used to discretize time and spatial terms in the given PDE. Since this contribution aims to propose a numerical method, its application to the flow problem is also given. Because of this study, more computational schemes can be constructed that will have more advantages than these schemes.

## 2 Numerical scheme

The proposed numerical scheme is a two-stage scheme consisting of explicit and implicit stages. The explicit stage is the forward Euler method, whereas the implicit stage comprises some unknown parameters. The first stage finds the problem’s solution at some unknown time level. The second stage of the scheme uses the information of the solution computed from the first stage of the scheme and finds a solution at

subject to the boundary conditions

where

where

The first stage of the scheme is given as,

This stage 4) requires only the information of the dependent variable computed at the previous time level. The next stage consists of a difference equation that contains some unknown parameters. The second stage of the scheme with unknown constant parameters

The unknown parameters

Substituting 4) and 6) into Eq. 5, it is obtained

Comparing coefficients of

By solving Eqs. 8–10, the values for

Therefore, the second stage of the proposed scheme for time discretization of Eq. 1 can be expressed as:

In Eq. 12,

Since the suggested scheme is only a time discretization scheme that can be used to discretize time derivative terms in considered time-dependent PDEs, it is to be noted that a semi-discrete scheme (12) can only be used to discretize time derivative terms. Space discretization can be performed by employing any scheme that discretizes spatial term(s) in given partial differential equations. For this contribution, a compact fourth-order scheme is adopted that can be used to discretize space derivative term(s). The adopted compact scheme requires finding the additional partial derivative of the given differential equation.

Applying a compact scheme to Eq. 1 using

where

Since Eq. 13 requires finding a fourth-order spatial derivative term, which will be found from Eq. 1 using

Substituting Eq. 14 into Eq. 13 gives

Equation 15 can be re-written as

Employing the first stage of the proposed scheme in Eq. 16 yields

Applying the second stage of the proposed scheme to Eq. 16, it follows

Since Eq. 18 contains second order time derivative term, that can be found by finding the second partial derivative of Eq. 1 using

Substituting Eq. 19 into Eq. 18, it follows

## 3 Stability analysis

The stability analysis for linear parabolic PDEs will be presented. This stability analysis gives the condition on the ratio of temporal step size and squared spatial step size, called diffusion number. So, a scheme with a large stability region might be preferred more than a smaller one. The scheme with a large stability region gives a large choice of step sizes and larger intervals for numerical values of involved parameters. So, the time consumed by the code for those schemes with large stability regions will be smaller than others. To start the stability analysis of the proposed scheme for classical or standard parabolic Eq. 1 using

Substituting some of the transformations from Eq. 21 into Eq. 17 yields

where

Dividing both sides of Eq. 22 by

Re-write Eq. 23 as

Substituting transformation from (21) into Eq. 20, it follows:

Dividing both sides of Eq. 25 by

Equation 26 can be re-written as

where

Since the proposed scheme is constructed on three-time levels, finding an amplification matrix requires one more equation. Also, note that the stability condition can be found without constructing an additional equation. But for this stability analysis, an additional condition is required. The additional equation is constructed as:

The vector matrix equation can be obtained from Eqs. 27, 28 as

The stability condition will be imposed on the Eigenvalues of the coefficients matrix of Eq. 29, which can be expressed as:

The stability conditions can be written as

So, if a scheme satisfies condition (31), it will be stable. A more detailed explanation of the stability condition can be obtained from a graphical view of condition (31). Those values of diffusion number

## 4 Problem formulation

This contribution also consisted of the mathematical model of boundary layer flow over the oscillatory sheet. The mathematical model will be comprised of a system of parabolic PDEs. One must first construct a matrix equation and discretize it using the compact scheme in the spatial direction and the proposed scheme in the temporal direction before determining the convergence condition for a system of PDEs. For this reason, consider a matrix vector equation of the form

where

Applying the first stage of the proposed scheme to Eq. 33, is obtained

Now applying the second stage of the proposed scheme on Eq. 33, yields

**Theorem 4.1:. **The scheme given in the form of Eqs. (34)-(35) converges.

**Proof:. **Let the first stage of the exact scheme be expressed as:

Let the error be

Re-write Eq. 37 as

Applying the norm on both sides of Eq. 38 and re-write the resulting equation in the form of

where

Subtracting Eq. 40 from Eq. 35 gives

Re-write Eq. 41, as

Applying norm on both sides of Eq. 42 and collecting the coefficients of

where

Let

where

Since the initial condition is exact so

Let

Let the error be bounded at first-time level

Substituting

Substituting

If this is continued, then for finite

Or

For infinite

## 5 Numerical examples

Two examples will be presented to check the scheme’s accuracy and convergence. The first example will be scalar linear parabolic PDE, and the second will be the system of parabolic equations arising in the boundary layer flow over the oscillatory sheet.

**Example 5.1:. ***Consider the following parabolic equation*

subject to the boundary conditions

and the initial condition is given as

The exact solution to the problem (53)-(55) is

**FIGURE 1**. Comparison of the proposed scheme with the Crank-Nicolson scheme for the norm of error over iterations using

**Example 5.2:. ***Casson NanoFluid Flow Over Oscillatory Sheet*Consider laminar, incompressible, unsteady, and unidirectional mixed convective Casson fluid flow over the oscillatory sheet. Let the sheet be moving with velocity

subject to the boundary conditions,

For reducing Eqs. 55–58 to dimensionless PDEs, the following transformation are considered:

Under transformations (60), Eqs. 56–59 can be expressed as:

subject to the dimensionless boundary conditions

The parameters are defined as follows:

For solving Eqs. 61–64, the proposed scheme in time direction with fourth-order compact scheme discretization in space is employed. But, the compact scheme requires the computation of a fourth-order spatial partial derivative, which is given as

Employing the fourth-order scheme to Eq. 61, which yields

Since the fourth-order spatial derivative is found in Eq. 65, so substituting Eq. 65 into Eq. 66 and the following equation can be obtained

Employing the first stage of the proposed scheme to Eq. 67, which is an explicit scheme, that gives

Employing the second stage of the proposed scheme on Eq. 67, which is an implicit stage of the scheme, yields

where

First, employ a compact scheme on Eq. 62 for spatial discretization. The semi-discrete equation can be expressed as

where

For time discretization, applying the first stage of the proposed scheme to Eq. 70, which gives

Employing the second stage of the proposed scheme to Eq. 70 gives

where

For getting a semi-discrete equation, applying a compact scheme to Eq. 63, is obtained

To get a fully discrete equation, employing the first stage of the proposed scheme to Eq. 73 and gives

Applying the second implicit stage of the proposed scheme on Eq. 80 gives

where

System of Eqs. 61–64 is solved using the proposed scheme with third-order accuracy in time. The first stage of the scheme can be applied to solve a given time-dependent parabolic equation without using any other scheme. But the second stage of the scheme requires using any other scheme to find a solution at the first time level. Two types of numerical schemes have been proposed. One of these schemes does not require any other scheme to find a solution at the first time level, but the second scheme requires either some estimated solution or a solution found by any different scheme. Mainly, scheme with higher accuracy converges faster than those having low order accuracy. So, high-order schemes may have been preferred over low-order schemes for finding solutions to parabolic equations with some chosen error norm. Since the proposed schemes are implicit, an additional iterative scheme is chosen for solving difference equation(s) obtained by applying the proposed scheme to given PDE(s). The considered iterative scheme for this contribution is the Gauss-Seidel iterative method. This iterative scheme requires some initial guess to begin finding a solution. The first iteration of the solution is found by considering the initial guess, and later on, the solution is found by utilizing previous consecutive iterations. The iterative procedure will stop if desired criteria are met. The given stopping criteria are based on the maximum norms of solutions for each equation in the system computed at two consecutive iterations. In this contribution, an iterative scheme is considered to solve difference equations obtained by applying the proposed scheme to parabolic equations. This process finds the solution at each grid point and time level. The difference equations can be solved exactly by forming a metric vector equation. But, in this study, only the iterative procedure is considered. Stability is one of the most important factors for getting the converged solution. But, the constructed Matlab code also shows the convergence and divergence of the solution that depends on the numerical values of chosen parameters and values of step sizes in some particular interval.Figure 2 deliberates comparing the proposed scheme with the existing Crank-Nicolson scheme. The maximum of norms for the solution of each Eqs. 61–64 is computed. The norm for the difference between two solutions computed at two consecutive iterations is found, and a maximum of three norms over each iteration is displayed in Figure 2. From this Figure 2, it can be observed that the proposed scheme converges faster than the existing scheme. Since the proposed scheme gives third-order accuracy in time and fourth-order in space, whereas Crank-Nicolson provides second-order accuracy in time and fourth-order in space, it converges faster than the existing scheme. Figs. 3–9 shows the effect of different parameters on velocity, temperature, and concentration profiles obtained by applying the proposed numerical scheme on Eqs. 61–64.Figure 3 shows the effect of the Casson parameter and magnetic parameter on the velocity profile. It is seen that the velocity profile declines by incrementing Casson and magnetic parameters. This fall in velocity profile due to the Casson parameter is the consequence of the diffusion coefficient reduction affecting the flow’s velocity. The rising values of the magnetic parameter diminish the flow’s velocity because the generation of Lorentz force in electrically conducting fluid resists the flow’s velocity, and consequently velocity of the flow drops. The effect of the porosity parameter and inertia coefficient on the flow’s velocity is discussed in Figure 4. Velocity profile de-escalates by enhancing the porosity parameter and inertia coefficient. Since the fluid’s viscosity resists the flow’s velocity, increasing the porosity parameter increases the fluid’s viscosity and reduces the velocity profile. Also, the velocity of the flow de-escalates by enhancing the inertia coefficient; this is the case when the drag force escalates, resulting in a decline in the flow’s velocity. The impact of thermal and solutal mixed convection parameters on the velocity profile is portrayed in Figure 5. Velocity profile rises by augmenting thermal and solutal mixed convection parameters. The rise in the velocity profile results from the temperature and concentration difference growth between the wall and the region away from the plate. Since the rise of temperature and concentration differences produces faster heat and concentration gradients in the flow, for mixed convection problems, temperature and concentration gradients are the factors responsible for generating flow in the fluid. Therefore, the flow’s velocity develops by growing thermal and solutal mixed convection parameters. The effect of the Eckert number and radiation parameter on the temperature profile is discussed in Figure 6. The rising Eckert number and radiation parameter boost the temperature profile. Since the temperature profile is affected by the dissipation due to internal friction of particles in the fluid flow and by improving Eckert number, the temperature profile escalates. Also, due to the entrance of radiations in the flow, heat flux increases; therefore, the temperature profile is boosted. Figure 7 displays the variation of thermophoresis and Brownian motion parameters on the temperature profile. The temperature profile is escalated by uprising thermophoresis and Brownian motion parameters. The increment in the thermophoresis parameter results in the growth of thermophoresis force, due to which heavy particles of fluid come closer to the plate, and lighter particles tend to move from the plate to its surroundings. This procedure of shifting particles improves the temperature of the flow. The boost in the Brownian motion parameter improves the random movement of particles in the fluid; therefore, the temperature profile is boosted. The influence of the Schmidt number and reaction rate parameter on the concentration profile is depicted in Figure 8. Concentration profile is de-escalated by rising Schmidt number and reaction rate parameter values. The increment in the Schmidt number produces low mass diffusivity, resulting in decay in the concentration profile. The augmentation in the reaction rate parameter yields the breaking and forming of particles that slow down the concentration profile. Figure 9 deliberates the variation of thermophoresis and Brownian motion parameters on the concentration profile. The concentration profile escalates by raising the values of the thermophoresis parameter, whereas it de-escalates by increasing the Brownian motion parameter values. Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 show the contour plots of velocity, temperature, and concentration profiles using sine boundary conditions. The variation of the oscillatory boundary can be seen in contour plots of velocity and temperature profiles. Figure 16 shows the comparison of the solution obtained by the proposed scheme with those obtained by the exact solution for solving Stokes’ second problem. The obtained results can be verified by looking at this comparison given in Figure 16.

**FIGURE 2**. Comparison of the proposed scheme with the Crank-Nicolson scheme for finding the maximum norms of errors over iterations using

**FIGURE 3**. Effect of Casson and magnetic parameters on velocity profile using **(A)** *M* = 0.5 **(B)** *β* = 3.

**FIGURE 4**. Effect of porosity parameter and inertia coefficient on velocity profile using **(A)** *F _{r}* = 0.1

**(B)**

*λ*= 0.1.

**FIGURE 5**. Effect of thermal and solutal mixed convection parameters on velocity profile using **(A)** *λ*_{2} = 0.1 **(B)** *λ*_{1} = 0.1.

**FIGURE 6**. Effect of Eckert number and radiation parameters on temperature profile using **(A)** *R _{d}* = 0.1

**(B)**

*E*= 0.5.

_{C}**FIGURE 7**. Effect of thermophoretic and Brownian motion parameters on temperature profile using **(A)** *N _{b}* = 0.15

**(B)**

*N*= 0.1.

_{t}**FIGURE 8**. Effect of Schmidt number and reaction rate parameter on concentration profile using **(A)** *γ* = 0.9 **(B)** *S _{c}* = 1.

**FIGURE 9**. Effect of thermophoretic and Brownian motion parameters on concentration profile using **(A)** *N _{b}* = 0.1

**(B)**

*N*= 0.1.

_{t}**FIGURE 10**. Contour plot of velocity profile using

**FIGURE 11**. Contour plot of temperature profile using

**FIGURE 12**. Contour plot of concentration profile using

**FIGURE 13**. Contour plot of velocity profile using

**FIGURE 14**. Contour plot of temperature profile using

**FIGURE 15**. Contour plot of concentration profile using

## 6 Conclusion

This contribution dealt with numerical schemes for time discretization of time-dependent PDEs. Two schemes have been proposed with third-order accuracy in time. According to von Neumann’s stability analysis and convergence of obtained solution from Matlab code, the scheme was conditionally stable. Moreover, a mathematical model for boundary layer Casson nanofluid flow over an oscillatory sheet has been constructed, influencing joule heating and chemical reaction. The dimensionless form of the model was solved with existing and proposed schemes. An iterative scheme was also adopted to solve discretized or difference equations obtained by proposing a scheme on the considered system of ODEs. Following this research, different applications for the current approach may be developed [48–52]. The main concluding points can be expressed as.

1) The proposed scheme converged faster than the existing Crank-Nicolson scheme.

2) Results have been validated for the velocity profile of the boundary layer flow phenomenon.

3) The velocity profile was decayed by enhancing Casson and magnetic parameters, and growing porosity parameter and inertia coefficient values inclined it.

4) Increasing values of thermophoresis and Brownian motion parameters raised the temperature profile.

5) The concentration profile was escalated by incrementing the thermophoresis parameter.

## Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

## Author contributions

Conceptualization, methodology, and analysis, YN; funding acquisition, WS; investigation, YN; methodology, MM; project administration, KA; resources, KA; supervision, MA.; visualization, KA; writing—review and editing, MA.; proofreading and editing, MM All authors have read and agreed to the published version of the manuscript.

## 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.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

1. Tagawa T. Numerical analysis of magnetohydrodynamic flows. *Fluids* (2020) 5(1):23. doi:10.3390/fluids5010023

2. Miura H. Extended magnetohydrodynamic simulations of decaying, homogeneous, approximately-isotropic and incompressible turbulence. *Fluids* (2019) 4:46 Accessed March 11, 2019. doi:10.3390/fluids4010046

3. Nazeer M, Saleem S, Hussain F, Zia Z, Khalid K, Feroz N. Heat transmission in a magnetohydrodynamic multiphase flow induced by metachronal propulsion through porous media with thermal radiation. *Proc Inst Mech Eng E: J Process Mech Eng* (2022) 095440892210752. doi:10.1177/09544089221075299

4. Bhattacharyya K. Dual solutions in boundary layer stagnation-point flow and mass transfer with chemical reaction past a stretching/shrinking sheet. *Int Commun Heat Mass Transfer* (2011) 38:917–22. doi:10.1016/j.icheatmasstransfer.2011.04.020

5. Bhattacharyya K. Dual solutions in unsteady stagnation-point flow over a shrinking sheet. *Chin Phys Lett* (2011) 28:084702. doi:10.1088/0256-307x/28/8/084702

6. Bhattacharyya K. MHD stagnation-point flow of Casson fluid and heat transfer over a stretching sheet with thermal radiation. *J thermodynamics* (2013) 2013:1–9. doi:10.1155/2013/169674

7. Bhattacharyya K. Heat transfer analysis in unsteady boundary layer stagnation-point flow towards a shrinking/stretching sheet. *Ain Shams Eng J* (2013) 4:259–64. doi:10.1016/j.asej.2012.07.002

8. Bhattacharyya K, Mukhopadhyay S, Layek G. Slip effects on an unsteady boundary layer stagnation-point flow and heat transfer towards a stretching sheet. *Chin Phys Lett* (2011) 28:094702. doi:10.1088/0256-307x/28/9/094702

9. You X, Xu H, Liao S. On the nonsimilarity boundary-layer flows of second-order fluid over a stretching sheet. *J Appl Mech* (2010) 77:021002. doi:10.1115/1.3173764

10. Liao S-J. A general approach to get series solution of non-similarity boundary-layer flows. *Commun Nonlinear Sci Numer Simulation* (2009) 14:2144–59. doi:10.1016/j.cnsns.2008.06.013

11. Nakhchi ME, Nobari M, Tabrizi HB. Non-similarity thermal boundary layer flow over a stretching flat plate. *Chin Phys Lett* (2012) 29:104703. doi:10.1088/0256-307x/29/10/104703

12. Kousar N, Liao S. Series solution of non-similarity natural convection boundary-layer flows over permeable vertical surface. *Sci China Phys Mech Astron* (2010) 53:360–8. doi:10.1007/s11433-010-0124-z

13. Kousar N, Mahmood R. Series solution of non-similarity boundary-layer flow in porous medium. *Appl Math* (2013) 4. doi:10.4236/am.2013.48A018

14. Hunegnaw D, Kishan N. Unsteady MHD heat and mass transfer flow over stretching sheet in porous medium with variable properties considering viscous dissipation and chemical reaction. *Chem Sci Int J* (2014) 4:901–17. doi:10.9734/acsj/2014/11972

15. Chu Y-M, Nazir U, Sohail M, Selim MM, Lee J-R. Enhancement in thermal energy and solute particles using hybrid nanoparticles by engaging activation energy and chemical reaction over a parabolic surface via finite element approach. *Fractal and Fractional* (2021) 5:119. doi:10.3390/fractalfract5030119

16. Abbas S. Z., Nayak M. K., Mabood F., Dogonchi A. S., Chu Y. M., Khan W. A. Darcy Forchheimer electromagnetic stretched flow of carbon nanotubes over an inclined cylinder: Entropy optimization and quartic chemical reaction. *Math Methods Appl Sci* (2020) 1–23. doi:10.1002/mma.6956

17. Liu C, Khan MU, Ramzan M, Chu Y-M, Kadry S, Malik MY, et al. Nonlinear radiative Maxwell nanofluid flow in a Darcy–Forchheimer permeable media over a stretching cylinder with chemical reaction and bioconvection. *Scientific Rep* (2021) 11:9391. doi:10.1038/s41598-021-88947-5

18. Al-Zubaidi A., Nazeer M., Subia G. S., Hussain F., Saleem S., Ghafar M. M. Mathematical modeling and simulation of MHD electro-osmotic flow of Jeffrey fluid in convergent geometry. *Waves in Random and Complex Media* (2021) 1–17. doi:10.1080/17455030.2021.2000672

19. Nazeer M, Hussain F, Khan MI, El-Zahar ER, Chu Y-M, Malik M, et al. Retracted: Theoretical study of MHD electro-osmotically flow of third-grade fluid in micro channel. *Appl Math Comput* (2022) 420:126868. doi:10.1016/j.amc.2021.126868

20. Nazeer M, Hussain F, Shabbir L, Saleem A, Khan MI, Malik M, et al. A comparative study of MHD fluid-particle suspension induced by metachronal wave under the effects of lubricated walls. *Int J Mod Phys B* (2021) 35:2150204. doi:10.1142/s0217979221502040

21. Ali N, Nazeer M, Javed T, Razzaq M. Finite element analysis of bi-viscosity fluid enclosed in a triangular cavity under thermal and magnetic effects. *The Eur Phys J Plus* (2019) 134:2–20. doi:10.1140/epjp/i2019-12448-x

22. Kumar A, Sugunamma V, Sandeep N. Impact of nonlinear radiation on MHD non-aligned stagnation point flow of micropolar fluid over a convective surface. *J Non-Equilibrium Thermodynamics* (2018) 43:327–45. doi:10.1515/jnet-2018-0022

23. Pramanik S. Casson fluid flow and heat transfer past an exponentially porous stretching surface in presence of thermal radiation. *Ain Shams Eng J* (2014) 5:205–12. doi:10.1016/j.asej.2013.05.003

24. Khan MI, Waqas H, Khan SU, Imran M, Chu Y-M, Abbasi A, et al. Slip flow of micropolar nanofluid over a porous rotating disk with motile microorganisms, nonlinear thermal radiation and activation energy. *Int Commun Heat Mass Transfer* (2021) 122:105161. doi:10.1016/j.icheatmasstransfer.2021.105161

25. Qayyum S, Khan MI, Masood F, Chu YM, Kadry S, Nazeer M. Interpretation of entropy generation in Williamson fluid flow with nonlinear thermal radiation and first-order velocity slip. *Math Methods Appl Sci* (2021) 44:7756–65. doi:10.1002/mma.6735

26. Song Y-Q, Khan SA, Imran M, Waqas H, Khan SU, Khan MI, et al. Applications of modified Darcy law and nonlinear thermal radiation in bioconvection flow of micropolar nanofluid over an off centered rotating disk. *Alexandria Eng J* (2021) 60:4607–18. doi:10.1016/j.aej.2021.03.053

27. Waqas H, Khan SA, Khan SU, Khan MI, Kadry S, Chu Y-M. Falkner-Skan time-dependent bioconvrction flow of cross nanofluid with nonlinear thermal radiation, activation energy and melting process. *Int Commun Heat Mass Transfer* (2021) 120:105028. doi:10.1016/j.icheatmasstransfer.2020.105028

28. Alaidrous AA, Eid MR. 3-D electromagnetic radiative non-Newtonian nanofluid flow with Joule heating and higher-order reactions in porous materials. *Scientific Rep* (2020) 10(1):14513–9. doi:10.1038/s41598-020-71543-4

29. Hafeez A, Khan M, Ahmed A, Ahmed J. Features of Cattaneo-Christov double diffusion theory on the flow of non-Newtonian Oldroyd-B nanofluid with Joule heating. *Appl Nanoscience* (2021) 1–8:265–72. doi:10.1007/s13204-020-01600-x

30. Khan A, Shah Z, Alzahrani E, Islam S. Entropy generation and thermal analysis for rotary motion of hydromagnetic Casson nanofluid past a rotating cylinder with Joule heating effect. *Int Commun Heat Mass Transfer* (2020) 119:104979. doi:10.1016/j.icheatmasstransfer.2020.104979

31. Kumar A, Tripathi R, Singh R, &Chaurasiya VK. Simultaneous effects of nonlinear thermal radiation and Joule heating on the flow of Williamson nanofluid with entropy generation. *Physica A: Stat Mech its Appl* (2020) 551:123972. doi:10.1016/j.physa.2019.123972

32. Khashi’ie NS, Arifin NM, Pop I, Wahid NS. Flow and heat transfer of hybrid nanofluid over a permeable shrinking cylinder with joule heating: A comparative analysis. *Alexandria Eng J* (2020) 59(3):1787–98. doi:10.1016/j.aej.2020.04.048

33. Uddin I, Ullah I, Ali R, Khan I, Nisar KS. Numerical analysis of nonlinear mixed convective MHD chemically reacting flow of Prandtl–Eyring nanofluids in the presence of activation energy and Joule heating. *J Therm Anal Calorim* (2021) 145(2):495–505. doi:10.1007/s10973-020-09574-2

34. Kazemi MA, Javanmard M, Taheri MH, Askari N. Heat transfer investigation of the fourthgrade non-Newtonian MHD fluid flow in a plane duct considering the viscous dissipation, joule heating and forced convection on the walls. *SN Appl Sci* (2020) 2(10):1752–14. doi:10.1007/s42452-020-03567-4

35. Ur Rasheed H, Saleem S, Islam S, Khan Z, Khan W, Firdous H, et al. Effects of joule heating and viscous dissipation on magnetohydrodynamic boundary layer flow of jeffrey nanofluid over a vertically stretching cylinder. *Coatings* (2021) 11(3):353. doi:10.3390/coatings11030353

36. Eswaramoorthi S, Loganathan K, Faisal M, Botmart T, Shah NA. Analytical and numerical investigation of Darcy-Forchheimer flow of a nonlinear-radiative non-Newtonian fluid over a Riga plate with entropy optimization. *Ain Shams Eng J* (2022) 14:101887. doi:10.1016/j.asej.2022.101887

37. Rasool G, Shah NA, El-Zahar ER, Wakif A. Numerical investigation of EMHD nanofluid flows over a convectively heated riga pattern positioned horizontally in a Darcy-forchheimer porous medium: Application of passive control strategy and generalized transfer laws. *Waves in Random and Complex Media* (2022).

38. Rauf A, Shah NA, Mushtaq A, Botmart T. Heat transport and magnetohydrodynamic hybrid micropolar ferrofluid flow over a non-linearly stretching sheet. *AIMS Math* (2023) 8(1):1–20. doi:10.1080/17455030.2022.2074571

39. Wakif A, Shah NA. Hydrothermal and mass impacts of azimuthal and transverse components of Lorentz forces on reacting Von Kármán nanofluid flows considering zero mass flux and convective heating conditions. *Waves in Random and Complex Media* (2022) 1–22.

40. Shah NA, Asogwa KK, Mahsud Y, Lee SR, Kang S, Chung JD, et al. Effect of generalized thermal transport on MHD free convection flows of nanofluids: A generalized atangana-baleanu derivative model. *Case Stud Therm Eng* (2022) 40:102480. doi:10.1080/17455030.2022.2136413

41. Zubair Akbar Qureshi M, Faisal M, Raza Q, Ali B, Botmart T, Shah NA. Morphological nanolayer impact on hybrid nanofluids flow due to dispersion of polymer/CNT matrix nanocomposite material. *AIMS Math* (2023) 8(1):633–56. doi:10.3934/math.2023030

42. Ibrahim A. Eigenvalue approach on fractional order theory of thermoelastic diffusion problem for an infinite elastic medium with a spherical cavity. *Appl Math Model* (2015) 39(20):6196–206. doi:10.1016/j.apm.2015.01.065

43. Abbas IA, Zenkour AM. Dual-phase-lag model on thermoelastic interactions in a semi-infinite medium subjected to a ramp-type heating. *J Comput Theor Nanoscience* (2014) 11(34):642–5. doi:10.1166/jctn.2014.3407

44. Priyadharshini P, Archana MV, Ahammad NA, Raju CSK, Yook SJ, Shah NA. Gradient descent machine learning regression for MHD flow: Metallurgy process. *Int Commun Heat Mass Transfer* (2022) 138:106307. doi:10.1016/j.icheatmasstransfer.2022.106307

45. Sajjan K, Shah NA, Ahammad NA, Raju CSK, Kumar MD, Weera W. Nonlinear Boussinesq and Rosseland approximations on 3D flow in an interruption of Ternary nanoparticles with various shapes of densities and conductivity properties. *AIMS Math* (2022) 7(10):18416–49. doi:10.3934/math.20221014

46. Zhang KZ, Shah NA, Vieru D, El-Zahar ER. Memory effects on conjugate buoyant convective transport of nanofluids in annular geometry: A generalized Cattaneo law of thermal flux. *Int Commun Heat Mass Transfer* (2022) 135:106138. doi:10.1016/j.icheatmasstransfer.2022.106138

47. Prabhakar Reddy B, Makinde OD, Hugo A. A computational study on diffusion-thermo and rotation effects on heat generated mixed convection flow of MHD Casson fluid past an oscillating porous plate. *Int Commun Heat Mass Transfer* (2022) 138:138 106389. doi:10.1016/j.icheatmasstransfer.2022.106389

48. Nawaz Y, Shoaib Arif M, Abodayeh K. An explicit-implicit numerical scheme for time fractional boundary layer flows. *Int J Numer Methods Fluids* (2022) 94:920. Wiley. doi:10.1002/fld.5078

49. Shoaib Arif M, Abodayeh K, Nawaz Y. The modified finite element method for heat and mass transfer of unsteady reacting flow with mixed convection. *Front Media S.A.* (2022) 10:952787. doi:10.3389/fphy.2022.952787

50. Nawaz Y, Arif MS, Shatanawi W, Nazeer A. An explicit fourth-order compact numerical scheme for heat transfer of boundary layer flow. *Energies* (2021) 14(12):3396. doi:10.3390/en14123396

51. Shatanawi W, Raza A, Arif MS, Rafiq M, Bibi M, Mohsin M. Essential features preserving dynamics of stochastic Dengue model. *Comp Model Eng Sci* (2021) 126(1):201–15. doi:10.32604/cmes.2021.012111

52 Arif MS, Raza A, Rafiq M, Bibi M, Abbasi JN, Nazeer A, Javed U. Numerical simulations for stochastic computer virus propagation model. *Comput. Mater. Contin* (2020) 62:61–77. doi:10.32604/cmc.2020.08595

## Nomenclature

Keywords: numerical schemes for ODEs, stability, casson fluid flow, joule heating and numerical solution, convergence and extension

Citation: Nawaz Y, Arif MS, Abodayeh K and Mansoor M (2023) Finite difference schemes for MHD mixed convective Darcy–forchheimer flow of Non-Newtonian fluid over oscillatory sheet: A computational study. *Front. Phys.* 11:1072296. doi: 10.3389/fphy.2023.1072296

Received: 17 October 2022; Accepted: 09 January 2023;

Published: 02 February 2023.

Edited by:

Antonio F. Miguel, University of Evora, PortugalReviewed by:

Nehad Ali Shah, Sejong University, Republic of KoreaAatef Hobiny, King Abdulaziz University, Saudi Arabia

Copyright © 2023 Nawaz, Arif, Abodayeh and Mansoor. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Muhammad Shoaib Arif, marif@psu.edu.sa