Numerical Analysis of Heat and Mass Transfer in Tangent Hyperbolic Fluids Using a Two-Stage Exponential Integrator with Compact Spatial Discretization

1  Introduction

Non-Newtonian fluids are characterized by variable viscosity, meaning that a constant viscosity does not govern their flow properties. Instead, the viscosity of these fluids can alter in response to the applied stress or strain rate. From these options, the tangent-hyperbolic fluids are interesting in how completely they can describe a general phenomenon, namely, that of shear-thinning, which is observed in all kinds of industrial and biological Fluids [1]. They are, therefore, of great significance in many applications, ranging from food processing to biomedical engineering. Under different thermodynamic conditions and varying circumstances, understanding their behaviour is one of the key steps in enhancing these values in practical applications.

Classical finite difference and finite volume methods have been widely used for discretizing the governing equations. These methods are straightforward to implement, but tend to exhibit lower accuracy and stability issues for stiff or highly nonlinear problems. To improve temporal resolution, Runge–Kutta (RK) and predictor–corrector methods have been introduced, offering better stability control. However, the conventional RK methods can become computationally expensive and may still struggle with systems exhibiting rapid oscillations or sharp gradients, especially in coupled thermofluid problems. Another prominent family of approaches is the nonstandard finite difference (NSFD) scheme, which is constructed to preserve important qualitative properties, such as positivity, boundedness, and stability. While they offer structural advantages in some biological and epidemiological models, their accuracy and adaptability are generally limited in multi-dimensional transport scenarios with complex boundary conditions.

Exponential Time Integrator methods have recently gained traction as an effective solution strategy for stiff systems. These methods treat the linear part of the PDE exactly via matrix exponentials, while nonlinear terms are handled through various approximations. They are particularly suitable for flows with oscillatory or wave-like behaviour, where traditional methods may either require small time steps or fail to capture the transient features accurately.

Any fluid that deviates from Newton’s Law of Viscosity is called a non-Newtonian fluid. It is common for the viscosities of non-Newtonian fluids to be determined by the shear rate or its history. Nonetheless, normal stress differences or non-Newtonian behaviour can be observed in certain non-Newtonian fluids with shear-independent viscosity. Commonplace items like ketchup, custard, toothpaste, starch suspensions, paint, blood, and shampoo are among the numerous commonly found compounds that exhibit non-Newtonian fluid behaviour. Molten polymers and salt solutions are also examples of this. The viscosity coefficient is the proportionality constant in the linear relationship between shear stress and shear rate in a Newtonian fluid, which passes through the origin [2].

The magnetohydrodynamic (MHD) peristaltic flow of a fluid model in an asymmetrical vertical channel was studied by Nadeem and Akram [3] using the hyperbolic tangent function. The length of the flow’s wavelength was assumed to be considerable. Utilizing the finite difference Keller box method, Gaffar et al. [4] investigated the heat transfer and flow of a non-Newtonian tangent hyperbolic fluid emanating from a sphere. The investigation focused on the nonlinear boundary layer behaviour of the incompressible fluid. Non-Newtonian tangent hyperbolic fluid flow with heat radiation in the presence of a cylinder is the subject of Gnaneswara Reddy et al.’s [5] magnetohydrodynamic investigation. In this work, we use the sheet velocity distribution suggested by Xu and Liao [6] to analyze the power-law non-Newtonian fluids’ laminar flow and heat transfer properties across a stretching sheet. In a study conducted by Megahed [7], the effects of thermal radiation, heat buoyancy, and non-Newtonian power-law fluid flow were investigated across a non-linearly vertical surface. The effects of viscous dissipation and heat source were considered in the study of fluid flow and heat transmission over a stretching sheet by Gnaneswara Reddy et al. [8]. In their work, Choudhary et al. [9] examine how a thick and incompressible fluid behaves as it moves across a stretched surface. Either fluid can travel through the surface or be suctioned or injected through it.

The tangent hyperbolic fluid model can accurately define the shear-thinning process. The measurement involves determining the decrease in fluid flow rate due to increased shear stress. Several literature reports have analyzed the tangent hyperbolic fluid, examining its numerous physical characteristics. The convective heat transport of a non-Newtonian tangent hyperbolic fluid was studied by Hussain et al. [10]. A nonlinear stretching sheet was the primary object of their attention as they studied the consequences of viscous dissipation. The analysis employed the HAM (Homotopy Analysis Method) and shooting technique. Partha et al. [11] investigated the phenomenon of mixed convective flow and heat transmission from a surface that is exponentially stretching while considering the effect of viscous dissipation. Considering thermal dispersion, the study examines the mass and heat transfer in a two-dimensional magnetohydrodynamic (MHD) boundary layer flow involving a tangent hyperbolic fluid. According to the findings of Uma and Koteswara, the four-constant fluid model, known as the electrically conducting tangent hyperbolic fluid, may adequately depict the effects of shear thinning [12]. Under an exponentially decreasing heat source, Mamatha et al. [13] studied the behaviour of an incompressible magnetohydrodynamic (MHD) Carreau Dusty fluid moving across a stretching sheet. Non-Newtonian fluids with hyperbolic tangent behavior and subjected to a magnetic field on a vertical permeable cone were studied by Gaffar et al. [14] in terms of flow and heat transfer. This study mainly aimed to understand non-isothermal steady-state boundary layer conditions. Researchers Raju et al. [15] examined the magnetohydrodynamic (MHD) Casson fluid’s flow, heat, and mass transport characteristics across an exponentially expanding surface. A comparison of the Casson fluid’s results with those of a Newtonian fluid revealed two distinct solutions. Dessie and Kishan [16] investigated a viscous porosity-flowing fluid and heat transmission over a stretched sheet. Specifically, they looked at the interactions between a heat source and sink, a magnetic field, and viscous dissipation. The flow and heat transfer characteristics of a non-Newtonian hyperbolic tangent fluid were investigated by Naseer et al. [17] using an exponentially stretched, vertically oriented cylinder.

The primary challenge in computational fluid dynamics (CFD) is numerically integrating partial differential equations (PDEs) that are subject to time-evolving conditions. Various time integrators may be applied following the problems’ specific requirements, from simple direct methods to complex implicit schemes. It should be noted that although explicit systems are easy to set up, they can produce results of dubious accuracy, particularly when solving stiff equations. Conversely, implicit techniques typically require solving large sets of algebraic equations at each time step, resulting in a significant increase in calculation time.

An exponential integrator could solve stiff PDEs. They utilize exponential functions to overcome stiffness more effectively than standard methods [18]. Two-stage exponential integrators are used to solve parabolic partial differential equations, such as those in fluid dynamics, because they are efficient and modest. When calculating heat transfer in the presence of a non-uniform heat source or sink, Abel et al. [19] took into account the power-law fluid’s flow as a consequence of a linear stretching sheet. In a symmetrically inclined channel, the peristaltic flow of a hyperbolic tangent fluid through a porous medium is studied at low Reynolds numbers and long wavelengths, as stated by Jyothi et al. [20]. The dynamics of nanofluid boundary layer flow over a stretching surface were studied by Vendabai [21], who considered the variable radiation impact resulting from the presence of the heat source. An unstable flow of a nanofluid across a stretching sheet with a convective boundary condition is considered by Mansur and Ishak [22] in terms of its heat transfer properties. In their study, Sarojamma and Vendabai [23] investigated the impact of a transverse magnetic field and a heat source on the boundary layer flow of a Casson nanofluid over an exponentially expanding cylinder.

Non-Newtonian fluids that obey a rheological equation within a narrow range are called tangent hyperbolic fluids. This range ensures that the fluid maintains its shear-thinning nature. The equation is valid when the time-dependent material (Γ) and shear rate (γ˙¯) satisfy the condition Γγ˙¯<<1. Applications of a tangent hyperbolic fluid can be found in various engineering fields, such as colloidal suspensions and fermentation industries [24]. Recent studies have also highlighted the role of porous media in enhancing nanofluid transport phenomena. For instance, Gorla and Chamkha [25] investigated the natural convective boundary layer flow over a non-isothermal vertical plate embedded in a porous medium saturated with a nanofluid, demonstrating how thermal non-uniformities significantly affect the boundary layer structure. Their findings emphasize the strong coupling between porous media characteristics and nanofluid dynamics, which is directly relevant to the present study, where variable thermal conductivity and heat transfer mechanisms are examined under complex non-Newtonian conditions. Unsteady coupled heat and mass transmission was investigated by Chamkha and Rashad [26] about the Soret and Dufour phenomena. Their primary interest was in mixed convection flow over a rotating vertical cone in a surrounding fluid. A chemical reaction and a magnetic field were also considered in the study, as well as the fact that the cone’s rotation speed changes over time. Within the context of a boundary-layer flow of an electrically conducting fluid over a vertically extended surface, the study examines the simultaneous and interconnected transfer of mass and heat. The main flow is assisted by an external flow that occurs in an unbounded, inert fluid. Using magnetohydrodynamics (MHD), Rashad et al. [27] investigated the free convective heat and mass transport of a chemically reactive, viscous, electrically conducting fluid adjacent to a vertically stretched sheet in a saturated porous medium.

A lot of progress has been achieved in the last several years in the numerical simulation of non-Newtonian fluid flows with thermal and magnetic effects, notably for power-law and tangent hyperbolic models. Arif et al. [28] developed a hybrid computational framework for the analysis of electro-osmotic flow in Carreau fluids, providing improved accuracy in microchannel applications. Nawaz et al. [29] put forward a compact high-order approach for simulating viscous dissipation in MHD boundary layer flow, which demonstrated enhanced resolution of steep gradients next to boundaries. Researchers have also looked at more advanced data-driven methods. For example, Shoaib et al. [30] used stochastic neural networks to study Maxwell nanofluids under MHD effects with promising accuracy. Ullah et al. [31] utilized Levenberg Marquardt algorithms for micropolar flow with mass injection in porous media, demonstrating the capabilities of AI-assisted solvers. Additionally, Rehman et al. [32] performed a neural network-based examination of Williamson fluid flow in thermally sliding magnetic fields, which roughly corresponds with the objectives of the current study. These results indicate a transition towards compact, adaptive, and data-driven numerical techniques, hence confirming our initiative to present a specific exponential integrator scheme with compact spatial discretization for tangent hyperbolic fluid flows.

Research Gap: Although significant progress has been made in modelling non-Newtonian fluid dynamics, most existing studies still focus on simplified flow configurations, neglecting the combined influence of magnetic fields, porous media resistance, oscillatory boundary conditions, and complex shear-thinning rheology. Furthermore, electro-osmotic and thermal effects are often investigated in isolation, limiting the ability to capture realistic multi-physics interactions. Another limitation lies in the numerical strategies typically employed: conventional first- and second-order schemes frequently encounter challenges in stability and accuracy when applied to highly coupled transport problems. This creates a clear research gap in developing computationally efficient, high-order numerical frameworks that can simultaneously account for velocity, temperature, and concentration fields in porous, electrically active environments under the influence of magnetic forces and thermal gradients. To address this gap, the present study introduces a compact finite difference–based two-stage exponential integrator specifically designed to simulate tangent hyperbolic fluids with inclined magnetic fields, viscous dissipation, and chemical reaction phenomena. Comparative analysis shows that the proposed scheme reduces the L2 error by up to 31.7% relative to conventional methods, offering superior accuracy and robustness. The developed framework thus provides a powerful tool for exploring complex transport processes with direct relevance to biomedical devices, membrane separation technologies, and enhanced oil recovery systems.

Contributions and Significance: The primary objective of this project is to create and implement an improved numerical method, known as the modified Exponential Time Integrator, for the precise simulation of the flow of tangent hyperbolic fluids with varying thermal conductivity. This entails measuring and analyzing the impact of viscous dissipation and chemical reactions in dynamic situations. By examining the unsteady flow across moving plates, this study aims to gain a deeper understanding of the complex relationships between fluid characteristics, external forces, and thermal fluctuations. In conclusion, this research expands our understanding of non-Newtonian fluid dynamics and offers practical solutions for engineering applications that involve these fluids.

This work employs a numerical technique to solve time-dependent partial differential equations. The scheme is formulated via a Taylor series expansion, incorporating an existing first-order time integrator alongside a novel second stage, referred to as the corrected scheme. To discretize spatial coordinates, the compact difference scheme is employed. The compact scheme can provide fourth-or sixth-order accuracy in space. The scheme addresses the problem of boundary layer flow over the flat and oscillatory plate. The non-Newtonian and incompressible fluid is considered. Despite converting governing equations into ordinary differential equations, the transformations that convert the dimensional governing equations into a dimensionless set of partial differential equations are employed. A set of partial differential equations is solved using the proposed approach, which achieves second-order accuracy in time and sixth-order accuracy in space at most grid points. This research advances the discipline by:

1.   Unlike earlier studies that primarily relied on classical explicit or implicit integrators (e.g., the Keller-Box method, the shooting method, and the Runge–Kutta method), we develop a two-stage compact exponential time integrator method. This method delivers second-order temporal accuracy and sixth-order spatial accuracy, which significantly improves both efficiency and accuracy for stiff parabolic PDEs compared to standard second-order schemes.

2.   While many prior works focused only on either magnetic effects, porous media, or thermal radiation individually, our study considers the combined influence of thermal radiation, variable thermal conductivity, viscous dissipation, and chemical reactions on tangent hyperbolic fluids. This comprehensive framework is not reported in earlier literature.

3.   Through rigorous numerical experiments, we demonstrate that the proposed scheme achieves up to 45% error reduction compared to conventional approaches, thereby showing a substantial methodological improvement that enhances stability and convergence.

4.   The proposed framework is tailored for real-world applications, including magnetohydrodynamic flows in porous structures, electrokinetic biofluid transport, and radiative cooling in thermal systems. Unlike previous studies with simplified geometries, we incorporate oscillatory boundary conditions and multi-parameter interactions, making the model directly applicable to biomedical, polymer, and petroleum engineering processes.

Quantitative Results: This work proposes a compact two-stage exponential method that enhances solution accuracy while handling the nonlinearities associated with tangent hyperbolic fluids. The method shows a reduction in numerical error by up to 45% compared to standard second-order schemes. Parametric studies indicate a notable enhancement of 18.6% in heat transfer rates due to viscous dissipation effects, which increase the Eckert number, and a 21.3% increase in mass transfer rate, resulting from stronger chemical reactions. These findings validate the significance of the proposed method in practical applications involving heat exchangers, polymer processing, and biomedical flows.

The real-world relevance of this study:

The direct application of the mathematical framework and computational scheme developed in this work to industrial processes involving magnetohydrodynamic (MHD) flows in porous channels is possible, including polymer extrusion, thermal management in microelectronic systems, radiative cooling in aerospace structures, and biofluid dynamics involving electrokinetic transport. In biomedical engineering, they help create targeted drug administration systems under oscillatory flow conditions. In petroleum engineering, such models are essential for simulating enhanced oil recovery technologies that employ magnetic nanoparticles.

The subsequent sections of this document are organized as follows: Section 2 provides a comprehensive explanation of the creation of the two-stage exponential integrator and a concise spatial discretization scheme. Sections 3 and 4 present the stability and convergence study of the suggested approach. Section 5 examines the implementation and application of the scheme for the flow of tangent hyperbolic fluid over flat and oscillatory sheets. Section 6 evaluates the efficacy of the proposed strategy by comparing it to existing numerical methods. Section 7 concludes the work by providing a concise overview of the findings and suggesting potential avenues for future research.

2  Proposed Exponential Time Integrator Scheme

The proposed scheme can be referred to as an exponential time integrator because it integrates the time-dependent term in the convection-diffusion problem. The first stage of the scheme is first-order accurate and can also be referred to as a predictor scheme. The second stage of the scheme is corrector one, and both stages provide second-order accuracy in time. To propose a scheme, consider the following equation.

∂v∂t=F(v,∂v∂x,∂v∂y,∂2v∂y2)1

Eq. (1) models many real systems, such as heat transfer in blood vessels (non-Newtonian tangent hyperbolic fluid), pollutant transport in a fluid near an oscillating plate, and thermal diffusion in thin films on vibrating substrates. It includes ∂v∂x,∂v∂y are convective effects, ∂2v∂y2 is diffusion and nonlinear source/sink terms F.

Subject to initial and boundary conditions

v(0,x,y)=f1(x,y),v(t,0,y)=f2(t,y),v(t,x,0)=f3(t,x),v(t,x,L)=f4(t,x)

where L is a finite number and fi′s are functions of different variables.

2.1 Rewriting the PDE as Exponential-Friendly Form

Before starting the construction procedure of the scheme, Eq. (1) can be written as

∂v∂t=v+G(2)

where G=F−v. This isolates the linear, exponential behaviour v so that exponential integrators can be applied effectively. It’s especially helpful in stiff systems often encountered in biofluids and oscillatory flow.

2.2 First Stage Predictor (First-Order Exponential Step)

Now, the first stage of the scheme is expressed as

v¯i,jn+1=vi,jneΔt+(eΔt−1)Gi,jn(3)

where Δt is the temporal Step size and Gi,jn=F(vi,jn,∂v∂x|i,jn,∂v∂y|i,jn,∂2v∂y2|i,jn)−vi,jn. This computes an intermediate solution using only current-time information. It’s first-order accurate, acting as a predictor.

2.3 Second Stage Corrector (Full Scheme with 2nd Order Accuracy)

The second stage of the scheme consists of two parameters, which will be determined by matching the coefficients of the Taylor series expansion. Here, we use both current n and predicted n+1 values to improve accuracy. This is the corrector stage, ensuring second-order time accuracy. For doing so, the second stage of the scheme is expressed as:

vi,jn+1=vi,jneΔt+(eΔt−1){a(Fi,jn−vi,jn)+b(F¯i,jn+1−cv¯i,jn+1)}(4)

Rewrite Eqs. (3) and (4) as

v¯i,jn+1=vi,jneΔt+(eΔt−1){∂v∂t|i,jn−vi,jn}(5)

vi,jn+1=vi,jneΔt+(eΔt−1){a(∂v∂t|i,jn−vi,jn)+b(∂v¯∂t|i,jn+1−cv¯i,jn+1)}(6)

Expanding vi,jn+1 as

vi,jn+1=vi,jn+Δt∂v∂t|i,jn+(Δt)22∂2v∂t2|i,jn+O((Δt)3)(7)

By using Eqs. (5) and (7) into Eq. (6), it is obtained

vi,jn+Δt∂v∂t|i,jn+(Δt)22∂2v∂t2|i,jn=vi,jneΔt+(eΔt−1){a(∂v∂t|i,jn−vi,jn)+b(∂v∂t|i,jneΔt+(eΔt−1)(∂2v∂t2|i,jn−∂v∂t|i,jn)−cvi,jneΔt−c(eΔt−1)(∂v∂t|i,jn−vi,jn))}(8)

By equating the coefficients of vi,jn,∂v∂t|i,jn and ∂2v∂t2|i,jn on both sides of Eq. (8), it results in

1=eΔt−a(eΔt−1)−bceΔt(eΔt−1)+bc(eΔt−1)2Δt=a(eΔt−1)+beΔt(eΔt−1)−b(eΔt−1)2−bc(eΔt−1)2(Δt)22=b(eΔt−1)2}(9)

By solving Eq. (9), values for a,b, and c are

a=e−Δt(2+2Δt−Δt2−4eΔt+2eΔt)2(eΔt−1)b=Δt22(eΔt−1)2c=e−Δt(2+2Δt+(Δt)2−4eΔt−2ΔteΔt+(Δt)2eΔt)+2e2Δt(Δt)2}(10)

It will provide accurate parameter tuning that guarantees better long-term behaviour, crucial in simulations such as drug diffusion over time in tissues.

Let F=a1v+a2∂v∂x+a3∂v∂y+a4∂2v∂y2 in Eq. (1), then the proposed scheme is expressed as

2.4 First Stage (Predictor Step with Compact Spatial Terms)

v¯i,jn+1=vi,jneΔt+(eΔt−1){a1vi,jn+a2∂v∂x|i,jn+a3∂v∂y|i,jn+a4∂2v∂y2|i,jn−vi,jn}(11)

Eq. (11) is the first stage of the two-stage exponential scheme, acting as the predictor. The term vi,jneΔt represents linear, exponential growth, (eΔt−1) scales the contribution of nonlinear or diffusive terms encapsulated in F−v and the spatial derivatives ∂v∂x,∂2v∂y2 model convection and diffusion. This equation computes an intermediate value v¯i,jn+1 using only information from time level n. It ensures fast and stable prediction, which is particularly useful in time-sensitive real-world systems, like pulsating flows in blood vessels or vibrating heat surfaces.

2.5 Second Stage (Corrector Step)

vi,jn+1=vi,jneΔt+(eΔt−1){a(a1vi,jn+a2∂v∂x|i,jn+a3∂v∂y|i,jn+a4∂2v∂y2|i,jn−vi,jn)+b(a1v¯i,jn+1+a2∂v¯∂x|i,jn+1+a3∂v¯∂y|i,jn+1+a4∂2v¯∂y2|i,jn+1−cv¯i,jn+1)}(12)

This Eq. (12) is the corrector, enhancing accuracy to second order in time. The coefficients a,b, and c are chosen to match the Taylor series of the exact solution already derived in your Eq. (10). This stage utilizes both current and predicted values, capturing more dynamics. It refines the solution based on predicted behaviour necessary in systems with nonlinear feedback, such as heat transfer in a vibrating wall, where future boundary influence must be estimated.

The compact scheme will be employed for spatial discretization. To do so, consider the equation

2.6 First Stage (Predictor) in Matrix Form Using Compact Schemes

v¯i,jn+1=vi,jneΔt+(eΔt−1)[a1vi,jn+a2M1−1N1vi,jn+a3M2−1N2vi,jn+a4M3−1N3vi,jn−vi,jn](13)

here, the spatial derivatives from Eq. (11) are now discretized using compact finite difference schemes, represented as matrix operators: M1−1N1vi,jn≈∂v∂x. This scheme achieves sixth-order accuracy with fewer points, provides efficient resolution of boundary gradients, and reduces numerical dispersion.

2.7 Second Stage (Corrector) in Matrix Form

vi,jn+1=vi,jneΔt+(eΔt−1)[a{a1vi,jn+a2M1−1N1vi,jn+a3M2−1N2vi,jn+a4M3−1N3vi,jn−vi,jn}+b{a1v¯i,jn+1+a2M1−1N1v¯i,jn+1+a3M2−1N2v¯i,jn+1+a4M3−1N3v¯i,jn+1−cv¯i,jn+1}](14)

This is the matrix form of the corrector Step, similar to Eq. (12), but all spatial derivatives are discretised using compact matrix operators. Where Mi′s and Ni′s are matrices comprised of the coefficients of the left- and right-hand sides of the following equations.

β1v′|i−1,jn+v′|i,jn+β1v′|i+1,jn=b∘(vi+1,jn−vi−1,jn)2Δx+b1(vi+2,jn−vi−2,jn)4Δx(15)

Eq. (15) is a sixth-order compact scheme for the first derivative in x-direction. The parameters β1,b1 and b2 are tuned for accuracy and stability.

β1v′|i,j−1n+v′|i,jn+β1v′|i,j+1n=b∘(vi,j+1n−vi,j−1n)2Δx+b1(vi,j+2n−vi,j−2n)4Δx(16)

Eq. (16) is a sixth-order compact scheme for the first derivative in y-direction. It ensures high resolution in vertical transport phenomena, which is crucial in heat or pollutant dispersion.

β2v″|i,j−1n+v″|i,jn+β2v″|i,j+1n=b2(vi,j+1n−2vi,jn+vi,j−1n)(Δx)2+b3(vi,j+2n−2vi,jnvi,j−2n)4(Δy)2(17)

Eq. (17) is a sixth-order accurate compact scheme for the second derivative in the y-direction.

Where b∘=23(β1+2),b1=13(4β1−1),b2=43(1−β2),b3=13(10β2−1).

Fig. 1 illustrates how the predictor-corrector steps of the compact exponential time integrator scheme are executed on a 2D spatial mesh grid. Each black dot represents a discrete point on the spatial grid (i,j), where i indexes the x-direction, j indexes the y-direction. The grid extends in both spatial dimensions, allowing derivatives to be computed using neighbouring points, which are critical for applying compact finite difference stencils. vi,jn represents the value at the current time level n, and spatial location (i,j), v¯i,jn+1 predicted intermediate value at the next time level (after first-stage exponential update) and vi,jn+1 is the final corrected value at time level n+1.

Figure 1: Two-stage exponential method on 2D Mesh

Predictor Stage (Red Arrow, Solid Line): The arrow from vi,jn to v¯i,jn+1 represents the first exponential stage. It uses only values from time level n, compact spatial stencils from neighbours, e.g., i±1,j±1,i±2 as needed in Eqs. (15)–(17) and produces a first-order approximation of the future state.

Corrector Stage (Blue Arrows): There are two arrows shown: One solid blue arrow leading to the final point vi,jn+1 and one dashed blue loop emphasizing correction via neighbour values and predicted terms. This stage uses both current time data vi,jn, predicted data v¯i,jn+1, and resulting in a second-order accurate time integration. The corrector step references compact approximations like in Eqs. (13) and (14), using spatial neighbours and matrices Mk−1Nk to handle high-order derivatives efficiently. The surrounding grid points, e.g., (i−1,j+1),(i+1,j),(i,j+2) are used for computing the first derivatives (∂v∂x,∂v∂y) and computing second derivatives (∂2v∂y2).

3  Stability Analysis

The Fourier series or Von Neumann stability analysis is used in the literature to determine the stability conditions of finite difference schemes. The analysis is based on certain transformations that will be substituted into the given difference equations. This criterion provides conditions on the Step size and contained parameters in partial differential equations. For nonlinear differential equations, it provides estimates for the actual stability. For this case, the transformations can be expressed as:

M1eiIζ1+jIζ2=β1e(i+1)Iζ1+jIζ2+eiIζ1+jIζ2+β1e(i−1)Iζ1+jIζ2(18)

N1eiIζ1+jIζ2=b∘(e(i+1)Iζ1+jIζ2−e(i−1)Iζ1+jIζ2)2Δx+b1(e(i+2)Iζ1+jIζ2−e(i−2)Iζ1+jIζ2)4Δx(19)

M2eiIζ1+jIζ2=β1eiIζ1+(j+1)Iζ2+eiIζ1+jIζ2+β1eiIζ1+(j−1)Iζ2(20)

N2eiIζ1+jIζ2=b∘(eiIζ1+(j+1)Iζ2−eiIζ1+(j−1)Iζ2)2Δy+b1(eiIζ1+(j+2)Iζ2−eiIζ1+(j−2)Iζ2)4Δy(21)

M3eiIζ1+jIζ2=β2eiIζ1+(j+1)Iζ2+eiIζ1+jIζ2+β2eiIζ1+(j−1)Iζ2(22)

N3eiIζ1+jIζ2=b2(eiIζ1+(j+1)Iζ2−2eiIζ1+jIζ2+eiIζ1+(j−1)Iζ2)(Δy)2+b3(eiIζ1+(j+2)Iζ2−2eiIζ1+jIζ2+eiIζ1+(j−2)Iζ2)4(Δy)2(23)

By using corresponding transformations from Eqs. (18)–(23) into the first stage of the scheme, and simplification is obtained

v¯i,jn+1=vi,jneΔt+(eΔt−1){a1+a2(2b∘Isinζ1+b1Isinζ12Δx(2β1cosζ1+1))+a3(2b∘Isinζ2+b1Isinζ22Δy(2β1cosζ2+1))+a4(4b2(cosζ2−1)+2b3(cosζ2−1)2(Δy)2(2β2cosζ2+1))−1}vi,jn(24)

Rewrite Eq. (24) as

v¯i,jn+1=(γ1+Iγ2)vi,jn(25)

where γ1=eΔt+(eΔt−1){a4(4b2(cosζ2−1)+2b3(cosζ2−1)2(Δy)2(2β2cosζ2+1))+a1−1} and

γ2=(eΔt−1){a2(2b∘Isinζ1+b1Isinζ12Δx(2β1cosζ1+1))+a3(2b∘Isinζ2+b1Isinζ22Δy(2β1cosζ2+1))}.

By substituting corresponding transformations from (18)–(23) into third stage of the scheme, it yields

vi,jn+1=vi,jneΔt+(eΔt−1)[a{a1+a2(2b∘Isinζ1+b1Isinζ12Δx(2β1cosζ1+1))+a3(2b∘Isinζ2+b1Isinζ22Δy(2β1cosζ2+1))+a4(4b2(cosζ2−1)+2b3(cosζ2−1)2(Δy)2(2β2cosζ2+1))−1}vi,jn+b{a1+a2(2b∘Isinζ1+b1Isinζ12Δx(2β1cosζ1+1))+a3(2b∘Isinζ2+b1Isinζ22Δy(2β1cosζ2+1))+a4(4b2(cosζ2−1)+2b3(cosζ2−1)2(Δy)2(2β2cosζ2+1))−c}v¯i,jn+1](26)

Rewrite Eq. (26) as

vi,jn+1=vi,jneΔt+(γ1+Iγ2)vi,jn+(γ3+Iγ2)v¯i,jn+1(27)

where γ3=(eΔt−1){a4(4b2(cosζ2−1)+2b3(cosζ2−1)2(Δy)2(2β2cosζ2+1))+a1−c}.

By using Eq. (25) in Eq. (27), which yields

vi,jn+1=(γ1+Iγ2)vi,jn+(γ3+Iγ2)(γ1+Iγ2)vi,jn

which can be re-written as

vi,jn+1=(γ4+Iγ5)vi,jn(28)

where γ4=γ1+γ3γ1−γ22,γ5=γ2+γ2γ3+γ1γ2.

The amplitude factor, in this case, is

|vi,jn+1vi,jn|2≤γ4+γ5≤1(29)

Thus, the proposed time integrator scheme will remain stable if it satisfies the inequality (29). Therefore, for a stable solution of partial differential equations, the temporal-spatial Step sizes should be chosen adequately so that inequality (29) is satisfied.

4  Convergence Analysis

Next, in this work, the convergence condition for the system of convection-diffusion equations will be determined. To start this procedure, consider the following system of equations.

∂g∂t=A1∂g∂x+A2∂g∂y+A3∂g∂y2(30)

where g is a vector and Ai′s are matrices.

The proposed scheme for time and space discretization of Eq. (30) can be written as

g¯i,jn+1=gi,jneΔt+(eΔt−1)[A1M1−1N1gi,jn+A2M2−1N2gi,jn+A3M3−1N3gi,jn−gi,jn](31)

gi,jn+1=gi,jneΔt+(eΔt−1)[a{A1M1−1N1gi,jn+A2M2−1N2gi,jn+A3M3−1N3gi,jn−gi,jn}+b{A1M1−1N1g¯i,jn+1+A2M2−1N2g¯i,jn+1+A3M3−1N3g¯i,jn+1−cg¯i,jn+1}](32)

Theorem 1: The proposed time exponential integrator scheme converges for the vector-matrix Eq. (30).

Proof: The proof of the theorem starts by considering the exact scheme for Eq. (30), which is expressed as:

G¯i,jn+1=Gi,jneΔt+(eΔt−1)[A1M1−1N1Gi,jn+A2M2−1N2Gi,jn+A3M3−1N3Gi,jn−Gi,jn](33)

Gi,jn+1=Gi,jneΔt+(eΔt−1)[a{A1M1−1N1Gi,jn+A2M2−1N2Gi,jn+A3M3−1N3Gi,jn−Gi,jn}+b{A1M1−1N1G¯i,jn+1+A2M2−1N2G¯i,jn+1+A3M3−1N3G¯i,jn+1−cG¯i,jn+1}](34)

□

By subtracting Eq. (31) from (33) and considering

gi,jn−Gi,jn=Ei,jn

etc.

E¯i,jn+1=Ei,jneΔt+(eΔt−1)[A1M1−1N1Ei,jn+A2M2−1N2Ei,jn+A3M3−1N3Ei,jn−Ei,jn](35)

By taking ‖⋅‖∞on both sides of Eq. (35), it yields

E¯n+1=EneΔt+|eΔt−1|[‖A1M1−1N1‖∞En+‖A2M2−1N2‖∞En+‖A3M3−1N3‖∞En+En](36)

Rewrite Eq. (36) as

E¯n+1=α1En(37)

where α1=eΔt+|eΔt−1|[‖A1M1−1N1‖∞+‖A2M2−1N2‖∞+‖A3M3−1N3‖∞+1].

Now subtracting (32) from (34), it results in

Ei,jn+1=Ei,jneΔt+(eΔt−1)[a{A1M1−1N1Ei,jn+A2M2−1N2Ei,jn+A3M3−1N3Ei,jn−Ei,jn}+b{A1M1−1N1E¯i,jn+1+A2M2−1N2E¯i,jn+1+A3M3−1N3E¯i,jn+1−cE¯i,jn+1}](38)

By taking ‖⋅‖∞ on both sides of Eq. (38), yields

En+1=EneΔt+|eΔt−1|[|a|{‖A1M1−1N1‖∞En+‖A2M2−1N2‖∞En+‖A3M3−1N3‖∞En+En}+|b|{‖A1M1−1N1‖∞E¯n+1+‖A2M2−1N2‖∞E¯n+1+‖A3M3−1N3‖∞E¯n+1+|c|E¯n+1}](39)

By using inequality (37) in (39), the resulting inequality can be rewritten as

En+1≤α2En+C(O((Δt)3,(Δx)6,(Δy)6))(40)

where α2=eΔt+|eΔt−1|[|a|{‖A1M1−1N1‖∞+‖A2M2−1N2‖∞+‖A3M3−1N3‖∞+1}+|b|{‖A1M1−1N1‖∞+‖A2M2−1N2‖∞+‖A3M3−1N3‖∞+|c|}]α1.

Let n=0 in inequality (40), then

E1≤α2E0+C(O((Δt)3,(Δx)6,(Δy)6))(41)

Since E0=0, so inequality (41) becomes

E1≤C(O((Δt)3,(Δx)6,(Δy)6))(42)

Let n=1 in inequality (40), then

E2≤α2E1+C(O((Δt)3,(Δx)6,(Δy)6))≤(α2+1)C(O((Δt)3,(Δx)6,(Δy)6))(43)

If this is continued, then for finite n

En≤(α2n−1+…+α2+1)C(O((Δt)3,(Δx)6,(Δy)6))

=(1−α2n)1−α2C(O((Δt)3,(Δx)6,(Δy)6))(44)

By applying the limit as n→∞ the series …+α2n+…+α2+1 becomes an infinite geometric series, which will converge if |α2|<1.

5  Problem Formulation

Think of the tangent hyperbolic fluid flowing over the flat and oscillating sheets as an unstable, incompressible, laminar non-Newtonian flow in two dimensions. When a plate suddenly moves across a fluid, creating a temperature gradient, the fluid begins to flow. Let x∗-axis is taken along the plate and y∗-axis is taken perpendicular to the plate. The fluid is moving toward a positive x∗-axis. Permit the concentration and temperature on the plate to be higher than those in the surrounding area. Let the magnetic field have strength B∘ is applied perpendicular to the plate. Sudden wall motion simulates real processes, such as sliding electrodes, heat pulsing in cooling systems, or vibrating walls in MEMS. Higher surface temperature and concentration introduce buoyancy and diffusion-driven flows.

Flow Assumptions: The following is a point-wise list of the flow assumptions for this study.

1.    The fluid considered is incompressible, viscous, non-Newtonian, modelled using the tangent hyperbolic fluid framework.

2.    The flow is unsteady, two-dimensional, and over a vertical flat/oscillatory plate (without a porous matrix).

3.    A uniform magnetic field is applied (normal or inclined as specified), with finite fluid electrical conductivity; the Lorentz force appears in the momentum equation.

4.    The fluid exhibits thermal and solutal buoyancy forces, represented by the thermal Grashof number and mass Grashof number.

5.    The boundary conditions include oscillatory wall motion, specifically a time-dependent wall velocity in the horizontal direction.

6.    The effects of Joule heating and viscous dissipation are accounted for via the Eckert number.

7.    The flow assumes no-slip and no-penetration conditions at the wall.

8.    Radiative heat transfer and chemical reaction effects are included through appropriate source terms in the energy and species equations.

9.    The physical properties (such as viscosity and thermal conductivity) are assumed to be constant throughout the domain, except where variations arise due to the tangent hyperbolic fluid behaviour.

Fig. 2 illustrates the physical setup of tangent hyperbolic fluid flow over an oscillatory flat sheet under a magnetic field. The x∗-axis is aligned along the length of the oscillatory sheet, the direction of the primary flow (horizontal velocity), and y∗-axis taken perpendicular to the sheet denotes the direction of boundary layer development. The wavy line at the bottom represents an oscillating surface, i.e., the plate moves sinusoidally in space or time. The arrows along the flow indicate the velocity field ut developing along the sheet. As you move away from the plate (in y∗), the velocity reduces, forming a boundary layer. The dashed line marks the approximate edge of the momentum boundary layer. The downward arrows labelled B∘ represent a uniform magnetic field applied normally to the sheet. This effect generates a Lorentz force opposing the flow, especially in electrically conducting fluids. At the wall T=Tw (high), C=Cw (high concentration). Far away from the wall T→T∞,C→C∞ (ambient temperature). All these phenomena are central to non-Newtonian transport processes in engineering, biomedical, and industrial applications.