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.

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 n + 1 t h time level. For proposing a scheme, consider the time-dependent partial differential equation in the form of ∂ u ∂ t = f u , ∂ 2 u ∂ x 2 , (1) subject to the boundary conditions u t , − L = α 1 ,   u t , L = α 2 , (2) where α 1 and α 2 are constants which have the same dimensions as the dependent variable u , and the initial condition can be given as u 0 , x = g x , (3) where g x is the function that only depends on x .

The first stage of the scheme is given as, u ¯ i n + 1 = u i n + Δ t ∂ u ∂ t i n , (4)

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 a 1 ,   a 2 ,   a 3   &   a 4 can be expressed as: u i n + 1 = u i n + Δ t a 1 ∂ u ∂ t i n + 1 + a 2 ∂ u ∂ t i n + a 3 ∂ u ¯ ∂ t i n + a 4 Δ t 2 ∂ 2 u ∂ t 2 i n , (5)

The unknown parameters a 1 , a 2 , a 3 and a 4 can be found by using Taylor series expansions for u i n + 1 and ∂ u ∂ t i n , as follows u i n + 1 = u i n + Δ t ∂ u ∂ t i n + Δ t 2 2 ∂ 2 u ∂ t 2 i n + Δ t 3 6 ∂ 3 u ∂ t 3 i n + O Δ t 4 , ∂ u ∂ t i n + 1 = ∂ u ∂ t i n + Δ t ∂ 2 u ∂ t 2 i n + Δ t 2 2 ∂ 3 u ∂ t 3 i n + O Δ t 3 . (6)

Substituting 4) and 6) into Eq. 5, it is obtained u i n + Δ t ∂ u ∂ t i n + Δ t 2 2 ∂ 2 u ∂ t 2 i n + Δ t 3 6 ∂ 3 u ∂ t 3 i n = u i n + Δ t a 1 ∂ u ∂ t i n + a 1 Δ t ∂ 2 u ∂ t 2 i n + a 1 Δ t 2 2 ∂ 3 u ∂ t 3 i n + a 2 ∂ u ∂ t i n + a 3 ∂ u ∂ t i n + a 3 Δ t ∂ 2 u ∂ t 2 i n + Δ t 2 a 4 ∂ 2 u ∂ t 2 i n . (7)

Comparing coefficients of a 1 , a 2 , a 3 and a 4 from Eq. 7, it is obtained 1 = a 1 + a 2 + a 3 , (8) 1 2 = a 1 + a 3 + a 4 , (9) 1 6 = a 1 2 . (10) By solving Eqs. 8–10, the values for a 1 , a 2 , a 3 and a 4 can be obtained as: a 1 = 1 3 ,   a 2 = 1 2 + a 4 , a 3 = 1 6 − a 4 . (11)

Therefore, the second stage of the proposed scheme for time discretization of Eq. 1 can be expressed as: u i n + 1 = u i n + Δ t 1 3 ∂ u ∂ t i n + 1 + 1 2 + a 4 ∂ u ∂ t i n + 1 6 − a 4 ∂ u ¯ ∂ t i n + a 4 Δ t 2 ∂ 2 u ∂ t 2 i n . (12) In Eq. 12, a 4 is a free parameter, and the suitable value of this parameter can be chosen so that the stable solution can be obtained.

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 f = ∂ 2 u ∂ x 2 , it gives ∂ u ∂ t i = δ x 2 u i − Δ x 2 12 ∂ 4 u ∂ x 4 i . (13) where δ x 2 u i = u i + 1 − 2 u i + u i − 1 Δ x 2 and Δ x is the spatial step.

Since Eq. 13 requires finding a fourth-order spatial derivative term, which will be found from Eq. 1 using f = ∂ 2 u ∂ x 2 as ∂ 4 u ∂ x 4 = ∂ ∂ t   δ x 2 u i . (14)

Substituting Eq. 14 into Eq. 13 gives ∂ u ∂ t i = δ x 2 u i − Δ x 2 12 ∂ ∂ t   δ x 2 u i . (15)

Equation 15 can be re-written as ∂ ∂ t u i + Δ x 2 12 ∂ ∂ t   δ x 2 u i = δ x 2 u i . (16)

Employing the first stage of the proposed scheme in Eq. 16 yields u ¯ i n + 1 + Δ x 2 12 δ x 2 u ¯ i n + 1 = u i n + Δ x 2 12 δ x 2 u i n + Δ t δ x 2 u i n . (17)

Applying the second stage of the proposed scheme to Eq. 16, it follows u i n + 1 + Δ x 2 12 δ x 2 u i n + 1 = u i n + Δ x 2 12 δ x 2 u i n + Δ t 1 3 δ x 2 u i n + 1 + 1 2 + a 4 δ x 2 u i n + 1 6 − a 4 δ x 2 u ¯ i n + 1 + Δ t a 4 ∂ 2 u ∂ t 2 i n . (18)

Since Eq. 18 contains second order time derivative term, that can be found by finding the second partial derivative of Eq. 1 using f = ∂ 2 u ∂ x 2 . Alternatively, this second-order partial time derivative in Eq. 18 can be numerically approximated by the second-order central difference formula as: ∂ 2 u ∂ t 2 i n = u i n + 1 − 2 u i n + u i n − 1 Δ t 2 (19)

Substituting Eq. 19 into Eq. 18, it follows u i n + 1 + Δ x 2 12 δ x 2 u i n + 1 = u i n + Δ x 2 12 δ x 2 u i n + Δ t { 1 3 δ x 2 u i n + 1 + 1 2 + a 4 δ x 2 u i n + 1 6 − a 4 δ x 2 u ¯ i n + 1 } + a 4 u i n + 1 − 2 u i n + u i n − 1 (20)

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 f = ∂ 2 u ∂ x 2 . Following transformations will be applied: u ¯ i n + 1 = E n + 1 e i I ψ ,   u i ± 1 n = E n e i ± 1 I ψ u i n + 1 = E n + 1 e i I ψ , u i ± 1 n + 1 = E n + 1 e i ± 1 I ψ u i n − 1 = E n − 1 e i I ψ (21)

Substituting some of the transformations from Eq. 21 into Eq. 17 yields E ¯ n + 1 e i I ψ + Δ x 2 12 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ Δ x 2 E n + 1 = E n e i I ψ + Δ x 2 12 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ Δ x 2 E n + Δ t Δ x 2 e i + 1 I ψ − 2 e i I ψ + e i + 1 I ψ E n (22) where I = − 1 .

Dividing both sides of Eq. 22 by e i I ψ , it yields E ¯ n + 1 1 + 1 6 cos ⁡ ψ − 1 = E n + 1 6 cos ⁡ ψ − 1 E n + 2 d cos ⁡ ψ − 1 E n . (23)

Re-write Eq. 23 as E ¯ n + 1 = 5 6 + 1 6 cos ⁡ ψ + 2 d cos ⁡ ψ − 1 E n 5 6 + 1 6 cos ⁡ ψ (24)

Substituting transformation from (21) into Eq. 20, it follows: E n + 1 e i I ψ + 1 12 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ E n + 1 = E n e i I ψ + 1 12 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ E n + d 1 3 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ E n + 1 + 1 2 + a 4 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ E n + 1 6 − a 4 e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ E ¯ n + 1 + a 4 E n + 1 − 2 E n + E n − 1 (25)

Dividing both sides of Eq. 25 by e i I ψ and substituting Eq. 24 into the resulting equation, it is obtained E n + 1 1 + 1 6 cos ⁡ ψ − 1 − 2 3 d cos ⁡ ψ − 1 a 4 = 1 + 1 6 cos ⁡ ψ − 1 + 2 d 1 2 + a 4 cos ⁡ ψ − 1 + 1 6 − a 4 cos ⁡ ψ − 1 5 6 + 1 6 cos ⁡ ψ + 2 d cos ⁡ ψ − 1 5 6 + 1 6 cos ⁡ ψ − 2 a 4 E n + a 4 E n − 1 . (26)

Equation 26 can be re-written as E n + 1 = a ¯ E n + b ¯ E n − 1 . (27) where a ¯ = 5 6 + 1 6 cos ⁡ ψ + 2 d 1 2 + a 4 cos ⁡ ψ − 1 + 1 6 − a 4 cos ⁡ ψ − 1 5 6 + cos ⁡ ψ 6 + 2 d cos ⁡ ψ − 1 5 6 + cos ⁡ ψ 6 − 2 a 4 5 6 + cos ⁡ ψ 6 − 2 3 d cos ⁡ ψ − 1 − a 4 , b ¯ = a 4 5 6 + cos ⁡ ψ 6 − 2 3 d cos ⁡ ψ − 1 − a 4

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: E n = E n + 0 E n − 1 . (28)

The vector matrix equation can be obtained from Eqs. 27, 28 as E n + 1 E n = a ¯ b ¯ 1 0 E n E n − 1 (29)

The stability condition will be imposed on the Eigenvalues of the coefficients matrix of Eq. 29, which can be expressed as: λ 1 = a ¯ − a ¯ 2 + 4 b ¯ 2 and   λ 2 = a ¯ + a ¯ 2 + 4 b ¯ 2 . (30)

The stability conditions can be written as λ 1 ≤ 1 and   λ 2 ≤ 1 . (31)

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 ″ d ″ that satisfy condition (31) can be chosen for converged solutions.

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 ∂ v ∂ t = A ∂ 2 v ∂ y 2 + B   ∂ v ∂ y + C v (32) where A ,   B , and C are 3 × 3 matrices, and v is a vector having three entries. By employing a fourth-order compact discretization to Eq. 31, the following semi-discretize scheme can be obtained: ∂ ∂ t v i + Δ y 2 12 δ y 2 v i + B A − 1 δ y v i = A + C Δ y 2 12 + Δ y 2 12 B A − 1 B δ y 2 v i + B + Δ y 2 12 B A − 1 C δ y v i + C v i . (33)

Applying the first stage of the proposed scheme to Eq. 33, is obtained v ¯ i n + 1 + Δ y 2 12 ( δ y 2 v ¯ i n + 1 + B A − 1 δ y v ¯ i n + 1 ) − v i n − Δ y 2 12 δ y 2 v i n + B A − 1 δ y v i n = Δ t A + C Δ y 2 12 + Δ y 2 12 B A − 1 B δ y 2 v i n + Δ t B + Δ y 2 12 B A − 1 C δ y v i n + Δ t C v i n (34)

Now applying the second stage of the proposed scheme on Eq. 33, yields v i n + 1 + Δ y 2 12 δ y 2 v i n + 1 + B A − 1 δ y v i n + 1 − v i n − Δ y 2 12 δ y 2 v i n + B A − 1 δ y v i n = Δ t a 1 A 1 δ y 2 v i n + 1 + B 1 δ y v i n + 1 + C v i n + 1 + b 1 A 1 δ y 2 v i n + B 1 δ y v i n + C v i n + c 1 A 1 δ y 2 v ¯ i n + 1 + B 1 δ y v ¯ i n + 1 + C v ¯ i n + 1 + e 1 v i n + 1 − 2 v i n + v i n − 1 (35)

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: V ¯ i n + 1 + Δ y 2 12 ( δ y 2 V ¯ i n + 1 + B A − 1 δ y V ¯ i n + 1 ) − V i n − Δ y 2 12 δ y 2 V i n + B A − 1 V i n = Δ t A 1 δ y 2 V i n + Δ t B 1 δ y V i n + Δ t C V i n . (36)

Let the error be v ¯ i n + 1 − V ¯ i n + 1 = e ¯ i n + 1 and subtracting Eq. 36 from Eq. 34, it is obtained e ¯ i n + 1 + Δ 2 12 ( δ y 2 e ¯ i n + 1 + B A − 1 δ y e ¯ i n + 1 ) − e i n − Δ y 2 12 δ y e i n + B A − 1 δ y e i n = Δ t A 1 δ y 2 e i n + Δ t B 1 δ y e i n + Δ t C e i n . (37)

Re-write Eq. 37 as 5 6 e ¯ n + 1 + Δ 2 12 ( e ¯ i + 1 n + 1 + e ¯ i − 1 n + 1 ( Δ y ) 2 + B A − 1 δ y e ¯ i n + 1 ) − e i n − Δ y 2 12 δ y 2 e i n + B A − 1 δ y e i n = Δ t A 1 δ y 2 e i n + Δ t B 1 δ y e i n + Δ t C e i n (38)

Applying the norm on both sides of Eq. 38 and re-write the resulting equation in the form of e ¯ n + 1 ≤ γ 1 e n , (39) where γ 1 = 1 + Δ y 2 12 ( 4 Δ y 2 + B A − 1 Δ y ) + 4 d A 1 + c B 1 + Δ t C 2 3 − Δ y 12 B A − 1 .

Now consider the second stage of the exact scheme for Eq. 32 can be expressed as: V i n + 1 + Δ y 2 12 δ y 2 V i n + 1 + B A − 1 δ y V i n + 1 − V i n − Δ y 2 12 δ y 2 V i n + B A − 1 δ y V i n = Δ t a 1 A 1 δ y 2 V i n + 1 + B 1 δ y V i n + 1 + C V i n + 1 + b 1 A 1 δ y 2 V i n + B 1 δ y V i n + C V i n + c 1 A 1 δ y 2 V ¯ i n + 1 + B 1 δ y V ¯ i n + 1 + C V ¯ i n + 1 + e 1 V i n + 1 − 2 V i n + V i n − 1 . (40)

Subtracting Eq. 40 from Eq. 35 gives e i n + 1 + Δ 2 12 δ y 2 e i n + 1 + B A − 1 δ y e i n + 1 − e i n − Δ y 2 12 δ y 2 e i n + B A − 1 δ y e i n = Δ t a 1 A 1 δ y 2 e i n + 1 + B 1 δ y e i n + 1 + C e i n + 1 + a 2 A 1 δ y 2 e i n + B 1 δ y e i n + C e i n + a 3 A 1 δ y 2 e ¯ i n + 1 + B 1 δ y e ¯ i n + 1 + C e ¯ i n + 1 + a 4 e i n + 1 − 2 e i n + e i n − 1 (41)

Re-write Eq. 41, as 5 6 I . D − a 1 d A 1 − a 4 e i n + 1 = − Δ y 2 12 e i + 1 n + 1 + e i − 1 n + 1 Δ y 2 + B A − 1 δ y e i n + 1 − e i n − Δ y 2 12 δ y 2 e i n + B A − 1 δ y e i n = Δ t a 1 A 1 δ y 2 e i n + 1 + B 1 δ y e i n + 1 + C e i n + 1 + a 2 A 1 δ y 2 e i n + B 1 δ y e i n + C e i n + a 3 A 1 δ y 2 e ¯ i n + 1 + B 1 δ y e ¯ i n + 1 + C e ¯ i n + 1 + a 4 − 2 e i n + e i n − 1 (42) Applying norm on both sides of Eq. 42 and collecting the coefficients of e i n + 1 on the left-hand side of the resulting equation, it yields e n + 1 ≤ γ 5 e n + γ 6 e n − 1 + M O Δ t 3 ,   Δ y 4 (43) where γ 5 = γ 3 γ 2 ,   γ 6 = γ 4 γ 2 and γ 2 = 5 6 I . D − a 1 d A 1 − a 4 − 1 3 − Δ y 12 B A − 1 − a 3 d a 1 A 1 − c a 1 B 1 γ 3 = 1 + Δ y 2 12 4 Δ y 2 + B A − 1 Δ y + d 4 a 2 A 1 + 4 a 3 A 1 γ 1 + c a 2 B 1 + a 3 B 1 γ 1 + Δ t a 2 C + a 3 C γ + 2 a 4 ,   γ 4 = a 4 ,   c = ∆ t ∆ x Let max e n , e n − 1 = e n − 1

Re-write inequality (43) as e n + 1 ≤ γ ¯ e n − 1 + M O Δ t 3 ,   Δ y 4 . (44) where γ ¯ = γ 5 + γ 6 .

Let n = 1 in (44). It yields e 2 ≤ γ ¯ e 0 + M O Δ t 3 ,   Δ y 4 (45) Since the initial condition is exact so e 0 = 0 , and therefore inequality (45) becomes, e 2 ≤ M O Δ t 3 ,   Δ y 4 (46)

Let n = 2 in (44), and this yields e 3 ≤ γ ¯ e 1 + M O Δ t 3 ,   Δ y 4 (47)

Let the error be bounded at first-time level i . e . e 1 ≤ B ¯ and (47) becomes e 3 ≤ γ ¯ B ¯ + M O Δ t 3 ,   Δ y 4 (48)

Substituting n = 3 in (44) yields e 4 ≤ γ ¯ e 2 + M O Δ t 3 ,   Δ y 4 ≤ 1 + γ ¯ M O Δ t 3 ,   Δ y 4 (49)

Substituting n = 4 in (44) gives e 5 ≤ γ ¯ e 3 + M O Δ t 3 ,   Δ y 4 ≤ γ ¯ 2 + 1 + γ ¯ M O Δ t 3 ,   Δ y 4 (50)

If this is continued, then for finite ″ m , ″ following one of the inequality is obtained e 2 m ≤ γ ¯ m − 1 + γ ¯ m − 2 + … + γ ¯ + 1 M O Δ t 3 ,   Δ y 4 = γ ¯ m − 1 γ ¯ − 1 M O Δ t 3 ,   Δ y 4 . (51) Or e 2 m ≤ γ ¯ m + γ ¯ m − 1 + γ ¯ m − 2 + … + γ ¯ + 1 M O Δ t 3 ,   Δ y 4 = γ ¯ m + γ ¯ m − 1 γ ¯ − 1 M O Δ t 3 ,   Δ y 4 . (52) For infinite ″ m ″ the infinite geometric series … + γ ¯ m + γ ¯ m − 1 + … + γ ¯ + 1 will converge if γ ¯ ≤ 1 .

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 ∂ u ∂ t = ∂ 2 u ∂ x 2 (53) subject to the boundary conditions u t , 0 = 0 = u t , π (54) and the initial condition is given as u 0 , x = sin x . (55)

The exact solution to the problem (53)-(55) is u t , x = e − t sin x . Equation 53 with boundary and initial conditions (54)–(55) is solved with the existing Crank-Nicolson scheme and proposed scheme. The space discretization is performed with a fourth-order compact scheme. The two schemes are compared in Figure 1. Figure 1 shows the error norm over iterations consumed by both schemes. Since the proposed scheme is more accurate than the existing Crank-Nicolson scheme, it converges faster than the existing scheme using the same step sizes for both schemes.

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 U w = U ∘ ⁡ cos ⁡ ⁡ a ω t * or U ∘ ⁡ sin ⁡ ⁡ a ω t * . Let the x * -axis be placed along the plate and y * -axis be placed perpendicular to the plate. Initially, the plate is at rest, and then it starts to move in positive x * -axis. Consider the effects of Joule heating, thermal radiation, and chemical reaction. Under the assumption of boundary layer theory, considering Darcy Forcheimmer’s characteristics and [47], the governing equations of the flow can be expressed as: ∂ u * ∂ t * = ν 1 + 1 β ∂ 2 u * ∂ y * 2 + g β ∘ T − T ∞ + β 1 T − T ∞ 2 + β 2 C − C ∞ + β 3 C − C ∞ 2 − σ B ∘ 2 ρ u * − ν k p u * − F u * 2 (56) ∂ T ∂ t * = α ∂ 2 T ∂ y * 2   + ν C p 1 + 1 β ∂ u * ∂ y * 2 + σ B ∘ 2 ρ C p u * 2 + τ D B ∂ C ∂ y * ∂ T ∂ y * + D T T ∞ ∂ T ∂ y * 2 (57) ∂ C ∂ t * = D B ∂ 2 C ∂ y * 2 + D T T ∞ ∂ 2 T ∂ y * 2 − k 1 C − C ∞ , (58) subject to the boundary conditions, u * t * ,   y * = U o cos a ω t *   o r   U o sin a ω t * ,   T t * ,   y * = T w ,   C t * , y * = C w   w h e n   y * = 0 , u * t * , y * → 0 ,   T T * , y * → 0 ,   C t * , y * → 0 ,   w h e n   y * → ∞ , (59)

For reducing Eqs. 55–58 to dimensionless PDEs, the following transformation are considered: u = u * U 0 , y = ω ν y * , t = ω t * , θ = T − T ∞ T w − T ∞ ,   ϕ = C − C ∞ C w − C ∞ (60)

Under transformations (60), Eqs. 56–59 can be expressed as: ∂ u ∂ t = 1 + 1 β ∂ 2 u ∂ y 2 − M u − K u − F r u 2 + λ 1 θ + λ 1 μ 1 θ 2 + λ 2 ϕ + λ 2 μ 2 ϕ 2 , (61) ∂ θ ∂ t = 1 P r 1 + 4 3 R d ∂ 2 θ ∂ y 2 + E c 1 + 1 β ∂ u ∂ y 2 + M E c u 2 + N b ∂ θ ∂ y ∂ ϕ ∂ y + N t ∂ θ ∂ y 2 (62) ∂ ϕ ∂ t = 1 S c ∂ 2 ϕ ∂ y 2 + 1 S c N t N b ∂ 2 θ ∂ y 2 − γ ϕ (63) subject to the dimensionless boundary conditions u t , y = cos a t o r ⁡ sin a t ,   θ t , y = 1 ,   ϕ t , y = 1   w h e n   y = 0 , u t , y → 0 ,   θ t , y → 0 , ϕ t , y → 0   w h e n   y → ∞ . (64)

The parameters are defined as follows: λ 1 = g β ∘ T w − T ∞ ω u ∘ ,   λ 2 = g β 2 C w − C ∞ ω u ∘ ,   M = σ B ∘ 2 ω ρ , K = ν ω k p , F r = c b ω k p 1 2 ,   P r = ν α , R d = 4 σ * T ∞ 3 k k *   ,   E c = U ∘ 2 c p T w − T ∞ ,   N b = τ D B C w − C ∞ ν , N t = τ D T T w − T ∞ ν T ∞ ,   S c = ν D B ,   γ = k 1 ω

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 ∂ 4 u ∂ y 4 = 1 1 + 1 β ∂ 3 u ∂ y 2 ∂ t + M ∂ 2 u ∂ y 2 + λ ∂ 2 u ∂ y 2 + F r u ¯ ∂ 2 u ∂ y 2 − λ 1 ∂ 2 θ ∂ y 2 − λ 1 μ 1 θ ¯ ∂ 2 θ ∂ y 2 − λ 2 ∂ 2 ϕ ∂ y 2   − λ 2 μ 2 ϕ ¯ ∂ 2 ϕ ∂ y 2 (65)