AI-Assisted Hybrid Solver for Skin Friction and Sherwood Number Prediction in Eyring–Prandtl Nanofluid Flow over a Riga Plate

1  Introduction

Non-Newtonian fluids play a central role in modern industrial, biomedical, and energy-related applications due to their complex rheological behaviour that cannot be modelled using classical Newtonian assumptions. Such fluids arise naturally in polymer melts, drilling muds, biofluids, paints, lubricants, and slurry transport systems. Among the many constitutive models proposed to describe non-Newtonian behaviour, the Eyring–Prandtl fluid model is recognized for its ability to capture shear-thinning effects arising from molecular activation-energy mechanisms. The Eyring-Prandtl formulation is an attractive choice for high-accuracy transport modelling, unlike power-law or Bingham models, because it can yield a physically realistic stress-strain relationship over a wide range of shear rates.

However, over the last few years, nanoparticles incorporated into base fluids to create nanofluids have attracted considerable interest due to their high heat and mass transfer rates. The nanoparticles increase thermal conductivity, diffusive transport, and chemical reactivity, which are essential in applications such as heat exchangers, cooling of electronic devices, chemical reactors, biomedical drug delivery, and solar energy systems. The joint play of non-Newtonian rheology and nanoparticle transport processes, specifically Brownian motion and Thermophoresis, results in highly non-linear transport behaviour, which requires sophisticated computational treatment.

One reason nanofluid flow is of interest is its interaction with electromagnetic actuation, where Lorentz forces can significantly affect the momentum boundary layers. In this regard, the Riga plate, initially proposed by Gailitis et al. [1], has proven to be a viable electromagnetic actuator for manipulating the boundary layer. The plate consists of periodically spaced permanent magnets and electrodes beneath the surface, which generate a Lorentz force parallel to the wall that can accelerate or decelerate the fluid layers near the wall, depending on the applied electrical conditions. The unmatched mechanism contributes to the plate of Riga, making it very efficient at reducing drag forces, stabilizing fluids in high-temperature systems, eliminating heat in high-temperature systems, micro-pumping, and improving mixing in microreactors. In turn, the study of the Eyring-Prandtl nanofluid flow over a Riga plate provides valuable insights into the underlying transport processes and design-oriented approaches for present-day electromechanical systems.

Choi [2] first used the phrase “nano liquid” to describe a solution containing nanoparticles suspended in a conventional fluid, such as water, oil, or ethylene glycol. In applications involving heat transfer, such as microprocessors, energy engines, hybrids, and temperature reduction, nanomaterials play an essential role. Two additional elements, Brownian motion and thermophoretic force, were added to Buongiorno’s [3] mathematical formulation. Then, via a thermal mechanism between two surfaces, Islam et al. [4] discussed the MHD flow of micropolar nanofluids. To understand how a nano liquid behaves, Alempour et al. [5] described an elliptical cross-section of tubes with spinning walls. Rasool et al. [6] analyzed the entropy of the nanoliquid’s Darcy-Forchheimer flow past the non-linear sheet. Tanveer and Malik [7] conducted numerical studies on the slip impact of a Ree-Eyring nanofluid in peristaltic flow across a curved surface. For instance, Khan et al. [8] investigated the influence of the Lorentz force on Ree–Eyring nanomaterial flow over a paraboloid stretching surface. They demonstrated that magnetic effects can effectively modulate velocity and thermal boundary layers. Nano-liquid flow has been extensively studied, with heat transfer and thermal conductivity as defining characteristics [9–13]. Mkhatshwa et al. [14] addressed the Hall effect and the effect of a chemical reaction on the physical properties of a hybrid nanofluid. An explanation of the three-dimensional Darcy-Forchheimer flow of an Eyring-Powell nanoliquid through a porous medium was provided by Tlili et al. [15]. Thermodynamic slip in Oldroyd-B fluid magnetized axi-symmetric flow with infinite stretchable disks was investigated by Khan et al. [16].

Over the last several decades, the inherent value of porous media has been recognized in an increasing number of important sectors. Porous media enhance nanoparticle scattering and increase the success rate by increasing the contact area between a fluid and a strong surface. Consequently, the effectiveness of conventional thermal frameworks is greatly enhanced [17]. When the flow rate is low (Reynolds number <1), Darcy’s law applies to flow through porous media. Even though Darcy’s law accurately describes most streams in porous media, it breaks down at high flow rates. Forchheimer correlations point to higher flow rates. The Forchheimer term, an additional non-linear velocity term, was first developed by Forchheimer [18].

At large Reynolds numbers, it was also described as a huge “Forchheimer term” by Muskat [19]. For example, Seddeek [20] looked at how viscous dissipation and thermophoresis affect a system, Francis et al. [21] looked at the Cattaneo-Christov heat transition model, Arulmozhi et al. [22] looked at how heat and mass transfer behave in MHD nanofluid flow over an infinite moving vertical plate under combined radiative and chemical reactive effects, Hayat et al. [23] looked at what happens when thermal conductivity changes, and every article that dealt with Darcy-Forchheimer flow was systematically searched. Using the Darcy-Forchheimer porous media model, Mohana and Kumar [24] successfully simulated three-dimensional hybrid nanofluid flow over a stretched sheet, incorporating form drag and convective heat transfer. In their study, Animasaun et al. [25] recorded the Darcy flow of a ternary-hybrid nanofluid. Prakash et al. [26] computed a Darcy-Forchheimer porous medium for an Electro-magneto-hydrodynamic (EMHD) ternary hybrid non-Newtonian nanofluid using zeta potential and wall suction/injection effects across a wedge.

Applying opposing electric and magnetic fields never results in a uniformly directed flow. A paradoxical electric current, the Hall current, flows through an electric field that is perpendicular to a magnetic field in the Hall effect [27]. One possible explanation is the Hall current, which is amplified by the strong magnetic field. When a fluid with electrical conductivity moves across a magnetic field, a phenomenon known as a transverse flow happens. This flow defies Ohm’s law. The Hall effect resolves several proven issues and makes it easier to identify various flow parameters in a flow field. Hall current’s effects on the maximum hydrodynamic drag (MHD) flow of a thick, non-compressible, electrically conducting fluid are of considerable theoretical and practical interest to scientists because of the many practical and theoretical uses for this phenomenon [28]. A few examples of the many applications of Hall current include flow-through filtering devices, thermal energy storage, Hall monitors, Hall accelerators, rotors, and MHD power sources [29]. Several devices rely on biomechanics and mechanics, such as biomicro pumps, bioreactors for membrane oxygenators, biochips, and magnetic resonance angiography (MRA) [30]. A nanofluid in a combined cone-and-disk geometry is analyzed using Hall currents, following the example by Basavarajappa and Bhatta [31]. Recent studies have explored heat and mass transfer in non-Newtonian Casson nanofluids over wedge-shaped stretching surfaces using advanced numerical techniques, highlighting the influence of nonlinear stretching and unsteady flow on transport characteristics [32]. Entropy generation and boundary layer behavior in Casson nanofluids have also been investigated under nonlinear and unsteady conditions, providing valuable benchmarks for modeling complex nanofluid systems [33].

Oil pipeline friction, central heating and refrigeration applications, enhanced capacity, and flow tracers are just a few examples of the many industrial, biological, chemical, and technological domains that make the transport and flow processes of non-Newtonian fluids fascinating to researchers and scientists [34]. Non-Newtonian materials have become more common in industrial sectors over the last few decades, while Newtonian materials have remained prevalent. Reason being, when compared to Newtonian fluids, non-Newtonian fluids are better at predicting elastic properties because of their higher viscoelastic qualities [35]. It is also widely recognized that most fluids in the body do not behave as Newtonian fluids. Some fluids that do not adhere to Newton’s rule of viscosity are actually viscoelastic due to their inherent properties [36]. It means that when shear tension is low, flow is absent. So, for fluids that aren’t Newtonian, yield stress is the sole stress. The list of non-Newtonian fluids includes things like concentrated juices, pulp, jellies, paints, sugar solutions, clay films, concentrated juices again, honey, drilling pozzolans, synthetic lubricants, speciality oils, and blood. A novel technology with enormous energy-saving potential, expandable consumption might be included in district heating and cooling systems, as well as building Heating, Ventilation, and Air Conditioning (HVAC) systems. Incorporating a colour scheme with non-Newtonian features yielded a highly effective flow tracker. To prevent dispersion and breakdown, they are used to create a tracer fluid. This fluid is subsequently introduced into a turbulent flow as a thin strip. The simple Navier-Stokes equations cannot sufficiently describe the flow field properties of non-Newtonian fluids because the mathematical formulation required to model the flow problem is complex. Through a penetrating exponentially vertical cone, Gomathi and Poulomi [37] explored the potential for optimizing the entropy of an EMHD Casson-Williamson penta-hybrid nanofluid.

A further level of complexity emerges when nanofluid flow is exposed to electromagnetic actuation, particularly in the presence of a Riga plate. Specialized electromagnetic surfaces include the Riga plate, which consists of alternating layers of permanent magnets and electrodes that generate a spatially varying Lorentz force parallel to the plate surface. This externally manageable force allows near-wall flow dynamics to be directly controlled, thereby enabling the effective manipulation of boundary layers, frictional resistance, and heat/mass transport rates. Riga plate configurations have been highly beneficial for drag reduction, flow stabilization, microfluidic pumping, thermal control at high temperatures, and nanocoating technologies. Therefore, the research on the Eyring-Prandtl nanofluid flow past a Riga plate is an essential contribution to both Theory and practical engineering of thermal mass transport.

The correct determination of the skin friction coefficient and the local Sherwood number is paramount for engineering. The skin friction coefficient measures surface drag and energy loss, and the Sherwood number measures the efficiency of mass transfer at the solid-fluid interface. Such amounts are fundamental to streamlining surface coating, separation technologies, biomedical mass diffusion, and catalytic surface reactions. However, under the combined effects of non-Newtonian rheology, MHD forces, nanoparticle diffusion, and electromagnetic control, the resulting governing equations become strongly non-linear and tightly coupled, rendering closed-form analytical solutions impractical.

Nanofluids are manufactured fluids made by dispersing nanoparticles into more traditional base fluids. These particles can be anything from metals and metal oxides to carbides and carbon nanotubes. Because of their improved thermophysical properties, these suspensions work wonders for increasing thermal conductivity. Nanofluids have many potential uses across industries, including hybrid power systems, home refrigeration, thermal management in electronics, medicinal procedures, and high-performance heat exchangers. Currently, three main coupling models are available. The first model, starting with the data-driven paradigm, allows the input-output link to be derived without any physical mechanisms. After that, the physical framework model uses AI algorithms to improve the previous models. The tertiary model combines physical processes with data-driven ones. While there are numerous AI algorithms available, data-driven models are where the ANN really shines. Aerodynamic modeling, turbulence modeling, specialized flows, and mass and heat transport phenomena are five areas where it has been effectively used. Physical model scenarios typically employ several other AI approaches, including recurrent neural networks (RNNs), naïve Bayes (NB) algorithms, and support vector machines (SVMs).

Recent studies indicate that machine learning models are being applied to estimate fluid thermophysical properties and to explore correlations between the skin friction, the Nusselt number, and other fluid mechanics variables. We look at some of the research. Shafiq et al. [38] first attempted to estimate the skin friction coefficient and Nusselt number using artificial neural networks (ANNs). In their analysis of the temperature-Biot number correlation, they found that ANN was the best solution. The Sherwood, Nusselt, and skin-friction coefficients have been approximated using an ANN model for non-Newtonian fluids over two spinning disks, as done by Zhao et al. [39]. A Darcy-Forchheimer mixed nanofluid flow in a revolving Riga disk was studied by Sharma et al. [40].

The prediction of hybrid nanofluid properties also uses an ANN in conjunction with the Levenberg-Marquardt backpropagation method. Predicting fluid properties was one of the ANN’s identified uses. Similarly, Boumari et al. [41] used ANNs and the Levenberg-Marquardt algorithm to investigate heat transfer through helical tubes and to calculate the Nusselt number. Two parameters were input to the model: the nanoparticle volume percentage and the Reynolds number. The exergy efficiency, frictional entropy, and heat entropy production of CoFe2O4/Water nanofluids were studied by Sundar and Mewada [42] using ANFIS and multilayer perceptrons as well. The impact of a spinning cylinder on the phase behavior variation in an L-shaped ventilation chamber was modelled by Ouri et al. [43] using an ANFIS-based model. When investigating the nanofluid Sherwood numbers in turbulent flows, Beiki [44] found that ANFIS was quite helpful, and when researching laminar flows, FCM-ANFIS was even more so.

Inspired by the structure and function of the real brain, artificial neural networks are computational representations of interconnected networks of neurons. To create forecasts using data patterns, explicit programming is superfluous. Predicting flows, simulating turbulence, estimating properties, enhancing systems, and controlling flows are all examples of fluid mechanics in action. Artificial neural networks (ANNs) aid these endeavours by developing turbulence models, accurately predicting flow behaviour, and optimizing system efficiency. Using ANN systems in conjunction with genetic algorithms, Ahmadi et al. [45] calculated pressure reductions in automotive radiators containing nanofluids. Using ANN modeling and multi-objective optimization, Hojjat [46] investigated nanofluids as coolants in shell-and-tube heat exchanger topologies. The Darcy-Forchheimer Ree-Eyring fluid flow with an activation energy and a convective boundary condition, as it pertains to a porous stretched layer, can be effectively simulated using an ANN, according to studies done by Shafiq et al. [47]. The constructed ANN is a practical engineering method that, as stated in [48], may deliver extremely accurate forecasts and can be used in reliability analysis via the inverse power law with Bayesian regularization. By incorporating nanofluid properties into their ANN-imitated microchannel heat exchanger model, Kamsuwan et al. [49] improved the simulation. To accurately anticipate heat transfer rates, fluctuations in the Nusselt number, and entropy generation, Habib et al. [50] created an ANN framework for heat and mass transfer through a porous Riga surface. The recent work [51], which uses a soft computing-based approach to enhance heat transfer in reactive nanofluid flows, aligns well with our focus on AI-assisted modeling of non-Newtonian flows with thermal mass transport mechanisms.

Moreover, simultaneous forecasting of the local Sherwood number and the skin friction coefficient is instrumental for describing the behaviour of surface drag and mass transfer. The parameters are fundamental for optimizing industrial equipment such as heat exchangers, nanocoating processes, separation equipment, biomedical sprays, and chemical reactors. Conventional analysis techniques fail to provide nonlinearities in coupled magnetohydrodynamic (MHD) forces, nanoparticle diffusion, and non-Newtonian rheology. Hence, to study such multi-physics interactions, numerical simulation is the best method of analysis. Yet numerical solvers are computationally intensive, particularly when parameter sweeps or optimization studies are required.

More recent studies have suggested moving beyond these constraints by adopting hybrid computational strategies that combine physics-based numerical workflows with machine learning (ML) and neural network (NN) models. Neural networks are beneficial for non-linear computations of flow parameters into engineering quantities of interest and can be trained quickly to provide rapid surrogate forecasts. NN-based models have been effectively used in fluid mechanics and heat transfer for predicting drag and thermal conductivity, and for predicting nanoparticle concentration and optimizing in complex geometries. Although these achievements have been made, their combination with non-linear and non-Newtonian nanofluid flows in the presence of electromagnetic surfaces, e.g., Riga plates, has been only briefly explored.

Consequently, high-accuracy numerical schemes have become essential for simulating such complex transport processes. Explicit and implicit Runge-Kutta methods, finite difference methods, and predictor-corrector methods, which are commonly used as classical time-integration methods, can be either unstable, inaccurate, or simply too expensive to solve stiff magnetohydrodynamic nanofluid systems. Exponential integrators are newer alternatives that have been effective in addressing these limitations in solving time-dependent partial differential equations because they are more stable and better at handling stiffness. However, even isolated exponential approaches can be inaccurate in solving non-linear source equations.

Both of these observations lead to the present work, which suggests a hybrid numerical scheme built based on an exponential integrator and a Runge-Kutta method. The obtained scheme is 3rd-order accurate, has three computational steps, and is designed to effectively handle first-order time-dependent non-linear partial differential equations governing the temporal dynamics of nanofluid transport. The stability and convergence behaviour of both scalar equations and coupled Partial differential equation (PDE) systems are rigorously analyzed, and the method is strongly mathematically reliable. It is followed by the application of the proposed scheme to simulate Eyring-Prandtl nanofluid flow over a Riga plate, accounting for coupled heat and mass transfer, thereby enabling the accurate calculation of velocity, concentration, the skin friction coefficient, and the Sherwood number profile.

Although the solutions to numerical simulations are high-fidelity, repeated simulations across a broad parameter space (Reynolds number, Prandtl number, Schmidt number, thermophoresis parameter, and Brownian motion parameter) may entail high computational costs. To address this challenge, the current study further incorporates artificial intelligence-based prediction using a neural network model. Recently, neural networks have been shown to possess an impressive ability to learn non-linear mappings between physical inputs and engineering outputs in the study of fluid mechanics, heat transfer, and transport modelling. Such models can make instant predictions with very little computational effort after training and are therefore ideal for real-time optimization and control.

Novelty of this study: In this work, a numerically generated dataset is used to train a neural network capable of predicting the skin friction coefficient and the local Sherwood number with high accuracy. To validate the neural model’s learning performance, error histograms and regression plots of predicted outputs vs. target data are presented, confirming the robustness and reliability of the AI-based predictor. This hybrid numerical–neural framework thus combines the physical interpretability of numerical simulation with the computational efficiency of artificial intelligence.

Finally, the performance of the proposed hybrid numerical scheme is assessed through a comparative analysis with an existing numerical method, demonstrating that the present scheme produces significantly lower numerical error for selected step sizes. This confirms both the efficiency and accuracy advantages of the developed approach. Here are the main points of this work’s contributions:

1.    We developed a third-order-accurate, three-stage hybrid exponential–Runge–Kutta numerical scheme for non-linear time-dependent PDEs.

2.    We provided a rigorous stability and convergence analysis for both scalar and coupled PDE systems.

3.    We developed mathematical models and numerical simulations of Eyring–Prandtl nanofluid flow over a Riga plate with heat and mass transfer.

4.    Accurate computation and AI-based prediction of skin friction coefficient and local Sherwood number.

5.    A comparative analysis demonstrates that the neural predictor achieves high accuracy while significantly reducing computational cost compared to repeated simulations.

6.    The combined framework serves as a promising hybrid modelling tool for real-time optimization of nanofluid-based thermal and mass transfer systems.

This study offers a novel collaboration between computational fluid dynamics and artificial intelligence to improve predictive capability and reduce computational burden in modelling complex non-Newtonian nanofluid flows. The results highlight the practical potential of hybrid models for advancing the design of next-generation engineering devices that involve MHD actuation and nanoscale transport processes.

In the remainder of this study, Section 2 describes the construction of the proposed hybrid numerical scheme, combining a modified exponential integrator with the classical Runge–Kutta method, and Section 3 details its convergence and stability. In Section 4, the Eyring-Prandtl nanofluid flow across a Riga plate is formulated in detail, accounting for electromagnetic forcing, Brownian motion, and thermophoresis on the velocity, temperature, and concentration profiles. The implementation of the scheme on benchmark problems and its validation against exact and MATLAB PDE solver solutions is discussed in Section 5. In Section 6, a neural network-based surrogate model is developed and trained using simulation data to predict the skin friction coefficient and local Sherwood number with high accuracy. Visualization of velocity, temperature, concentration distributions, error contours, and AI regression plots is presented in Section 6 to support the findings. Finally, Section 7 offers concluding remarks, summarising the key insights and highlighting the practical relevance of the proposed hybrid computational-intelligent framework for real-time prediction of non-Newtonian nanofluid transport phenomena.

2  Proposed Predictor Corrector Scheme

A numerical scheme is constructed for discretizing time-dependent partial differential equations. The scheme only discretizes time-dependent terms. The first two stages of the scheme are called predictor stages, while the last stage is called the corrector stage. Before applying a scheme to a differential equation domain, the domain is divided into a finite number of equal subintervals. The solution will be found at the point of each subinterval. We split the domain into a space grid and a time grid: Space: y∈[0,L], divided into Ny−1 equal subintervals, grid points yj, spacing Δy. Time: t∈[0,tf], divided into Nt−1 equal time steps, grid points tn, step size Δt. We will compute approximations: vjn≈v(tn,yj). The differential equation can be expressed as:

∂v∂t=α1∂v∂y+α2∂2v∂y2.(1)

Subject to the boundary conditions

v(t,0)=f1(t),v(t,L)=f2(t).(2)

The initial condition is expressed as

v(0,y)=f3(y),(3)

where f1,f2,f3 are functions, v(t,y) could be velocity, temperature, or concentration, depending on the physical problem, α1∂v∂y is a convection (advection) term: transport due to bulk motion (e.g., fluid moving along a plate) and α2∂2v∂y2 is a diffusion term: smoothing due to molecular/Brownian diffusion or viscous effects. So, we are essentially designing a scheme that can accurately track how velocity or concentration evolves over time and space under combined convection–diffusion effects.

2.1 Stage 1: Exponential-Type Predictor

The first Stage of the proposed scheme is expressed as:

v¯jn+1=12vjne2Δt+12vjne−2Δt+(e2Δt+e−Δt−2){∂v∂t|jn−12vjn}.(4)

This scheme is first-order accurate in time. This is an exponential integrator-type step. It uses information at time level n: vjn, ∂v∂t|jn (which itself depends on spatial derivatives). It’s first-order accurate in time and used as a cheap predictor. Stage 1 provides a fast, rough forecast of how the field (velocity, concentration) will evolve over the next time step, accounting for the current local time derivative.

Physical Significance of Eq. (4): Eq. (4) represents the first-stage exponential predictor in the proposed three-stage time integration scheme. This step provides a computationally inexpensive, explicit forecast of the solution at the next time level using the known values of the solution and its time derivative at the current time level n. It belongs to the family of exponential integrators, which are particularly useful in stiff PDEs (e.g., magneto-convective non-Newtonian flows) where conventional explicit methods may require impractically small time steps. The exponential terms in Eq. (4) are introduced to enhance stability and accuracy beyond that of a standard forward Euler predictor. This formulation approximates the temporal evolution in a nonlinear system more accurately when the PDE contains sharp gradients or strong source terms (e.g., electromagnetic force or viscous dissipation).

2.2 Stage 2: Improved Predictor

The second Stage of the scheme is given as:

v¯¯jn+1=avjn+bv¯jn+1+14Δt∂v¯∂t|jn+1.(5)

It combines the old solution vjn and the first predictor v¯jn+1 and uses a new time derivative ∂v¯∂t evaluated at time n+1 (but with the predictor). The coefficients a,b are not arbitrary; they will be chosen (with c,d) to enforce third-order accuracy via Taylor Expansion. Stage 2 refines the forecast using improved information about how the PDE responds at the new time level, which is important in systems where the dynamics are sensitive (e.g., strong shear, strong diffusion, electromagnetic forcing).

2.3 Stage 3: Corrector

The last Stage of the scheme is expressed as:

vjn+1=cvjn+dv¯¯jn+1+23Δt∂v¯¯∂t|jn+1.(6)

This step corrects the solution using the original state vjn, the improved predictor v¯¯jn+1, and the updated time derivative. Again, c,d are chosen so that the final scheme matches the Taylor series of the exact solution up to O(Δt3). For our nanofluid flow, it ensures that boundary-layer development, velocity gradients, and concentration gradients are accurately tracked over time. This accuracy is crucial when we compute skin friction and Sherwood number, which depend on derivatives at the wall. The unknown parameters a,b,c, and d will be determined by applying the Taylor series Expansion. For doing this, substitute the first stage (4) into the second stage of the scheme (5) as:

v¯¯jn+1=avjn+b[12vjne2Δt+12vjne−2Δt+(e2Δt+e−Δt−2){∂v∂t|jn−12vjn}]+14Δt[12e2Δt∂v∂t|jn+12e−2Δt∂v∂t|jn+(e2Δt+e−Δt−2){∂2v∂t2|jn−12∂v∂t|jn}].(7)

Eq. (7) can be expressed as

v¯¯jn+1=(a+b)vjn+b(e2Δt+e−Δt−2)∂v∂t|jn+14Δt{∂v∂t|jn+(e2Δt+e−Δt−2)∂2v∂t2|jn}.(8)

By putting Eq. (8) into Eq. (6), it yields

vjn+1=cvjn+d[(a+b)vjn+b(e2Δt+e−Δt−2)∂v∂t|jn+14Δt{∂v∂t|jn+(e2Δt+e−Δt−2)∂2v∂t2|jn}]+23Δt[(a+b)∂v∂t|jn+b(e2Δt+e−Δt−2)∂2v∂t2|jn+14Δt{∂2v∂t2|jn+(e2Δt+e−Δt−2)∂3v∂t3|jn}].(9)

Rewrite Eq. (9) as

vjn+1=[c+d(a+b)]vjn +[bd(e2Δt+e−Δt−2)+d4Δt+23Δt(a+b)] ∂v∂t|jn]+[dΔt4(e2Δt+e−Δt−2)+2b3Δt(e2Δt+e−Δt−2)+(Δt)26∂2v∂t2|jn]+[(Δt)26(e2Δt+e−Δt−2)∂3v∂t3|jn].(10)

Now expanding vjn+1 using the Taylor series as

vjn+1=vjn+Δt∂v∂t|jn+(Δt)22∂2v∂t2|jn+(Δt)36∂3v∂t3|jn+O((Δt)4).(11)

Using Taylor series Expansion (11) in Eq. (10) and comparing the coefficients of vjn,∂v∂t|jn,∂2v∂t2|jn and ∂3v∂t3|jn on both sides of the resulting equations, it gives

c+d(a+b)=1bd(e2Δt+e−Δt−2)+d4Δt+23Δt(a+b)=ΔtdΔt4(e2Δt+e−Δt−2)+2b3Δt(e2Δt+e−Δt−2)+(Δt)26=(Δt)22(Δt)26(e2Δt+e−Δt−2)=(Δt)36}.(12)

Solving a system of Eq. (12) yields

a=116(16−12d+9d2),b=18(4−3d),c=116(16−24d+18d2−9d3) and d is a free parameter. This matching guarantees that, for reasonably smooth solutions, the method is third-order accurate in time, meaning that the time integration error is very small even for moderately large Δt. This is particularly valuable for stiff problems such as non-Newtonian MHD flows, where simple explicit schemes require extremely small time steps to remain stable.

2.4 Compact Spatial Discretization

Up to now, we treated ∂v∂t as if we know it exactly. In practice,∂v∂t|jn=α1∂v∂t|jn+α2∂2v∂y2|jn. Since the proposed scheme discretizes only the time derivative in Eq. (1), a compact scheme is employed for space discretization. Therefore, the proposed scheme using compact spatial discretization is given as

v¯jn+1=12vjne2Δt+12vjne−Δt+(e2Δt+e−Δt−2){α1A1−1B1vjn+α2A2−1B2vjn−12vjn},(13)

v=jn+1=avjn+bv¯jn+1+14Δt{α1A1−1B1v¯jn+1+α2A2−1B2v¯jn+1},(14)

vjn+1=cvjn+dv¯¯jn+1+23Δt{α1A1−1B1v¯¯jn+1+α2A2−1B2v¯¯jn+1},(15)

where A1,A2 and B1,B2 are matrices that are constructed from the left-hand and right-hand sides of the following equations, respectively.

For first derivative:

β1v′|j+1n+v′|in+β1v′|j−1n=b∘(vj+1n−vj−1n)2Δy+b1(vj+2n−vj−2n)4Δy.(16)

For the second derivative:

β2v′′|j+1n+v′′|jn+β2v′′|j−1n=b2(vj+1n−2vjn+vj−1n)(Δy)2+b3(vj+2n−2vjn+vj−2n)4(Δy)2.(17)

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

These relations can be written in matrix form: A1vy=B1v, so vy≈A1−1B1v. A2vyy=B2v, so vyy≈A2−1B2v. So, in our final scheme: ∂v∂t|jn≈α1A1−1B1vn+α2A2−1B2vn.

This scheme gives a higher-order spatial accuracy using short stencils (involving only a few neighbouring points). It captures sharp gradients better, which is essential when we calculate wall shear stress (forskin friction) and the wall concentration gradient (for the Sherwood number).

3  Stability and Convergence Analysis

To determine the necessary conditions for the stability of finite difference schemes, authors have used either the von Neumann or the Fourier series stability analyses. The analysis transforms the difference equations into a trigonometric equation. The analysis considers both time- and space-discretization schemes. It gives exact stability conditions of finite difference schemes applied to a linear differential equation. It may apply to a non-linear differential equation after linearization. For non-linear differential equations, it provides an estimate of the exact stability conditions. The transformations for this case can be expressed as:

A1eiIψ=β1e(j+1)iψ+eiIψ+β1e(j−1)iψ,(18)

B1eiIψ=b∘(e(j+1)iψ−e(j−1)Iψ)2Δy+b1(e(j+2)iψ−e(j−2)iψ)4Δy,(19)

A2eiIψ=β2e(j+1)iψ+eiIψ+β2e(j−1)iψ,(20)

B2eiIψ=b2(e(j+1)iψ−2eiIψ+e(j−1)iψ)(Δy)2+b3(e(j+2)iψ−2eiIψ+e(j−2)iψ)4(Δy)2,(21)

where i=−1 is an imaginary number.

Using transformations (18)–(21) into the first stage of the proposed scheme (13) gives

v¯jn+1=12vine2Δt+12vine−Δt+(e2Δt+e−Δt−2){α1(2b∘isinψ+b1isin2ψ2Δy(2β1cos⁡ψ+1))+α2(4b2(cosψ−1)+b3(cos2ψ−1)2(Δy)2(2β2cos⁡ψ+1))−12}vjn.(22)

Eq. (22) can be re-written as

v¯jn+1=(c1+ic2)vjn,(23)

where c1=12e2Δt+12e−Δt+(e2Δt+e−Δt−2){α2(4b2(cosψ−1)+b3(cos2ψ−1)2(Δy)2(2β2cos⁡ψ+1))−12}.

And c2=(e2Δt+e−Δt−2){α1(2b∘sinψ+b1sin2ψ2Δy(2β1cos⁡ψ+1))}.

Now using transformations (18)–(21) into the second stage of the scheme (14) yields.

v=jn+1=avjn+bv¯jn+1+14Δt{α1(2b∘isinψ+b1isin2ψ2Δy(2β1cos⁡ψ+1))+α2(4b2(cosψ−1)+b3(cos2ψ−1)2(Δy)2(2β2cos⁡ψ+1))}v¯jn+1.(24)

Rewrite Eq. (24) as

v=jn+1=avjn+(c3+ic4)v¯jn+1,(25)

where c3=b+14Δtα1(2b∘sinψ+b1sin2ψ2Δy(2β1cos⁡ψ+1)) and c4=14Δtα2(4b2(cosψ−1)+b3(cos2ψ−1)2(Δy)2(2β2cos⁡ψ+1)).

Substituting Eq. (23) into Eq. (25) leads to

v=jn+1=(c5+ic6)vjn,(26)

where c5=a+c1c3−c2c4 and c6=c1c4+c2c3.

Now using transformations (18)–(21) into the third stage of the scheme (15) is

vjn+1=cvjn+(c7+ic8)v=jn+1.(27)

Substituting Eq. (26) into Eq. (27) yields

vjn+1=(c9+ic10)vjn,(28)

where c9=c+c5c7−c6c8 and c10=c5c8+c6c7,

where c7=d+2Δt3α1(2b∘sinψ+b1sin2ψ2Δy(2β1cos⁡ψ+1)) and c8=2Δt3α2(4b2(cosψ−1)+b3(cos2ψ−1)2(Δy)2(2β2cos⁡ψ+1))

The amplification factor for this case is written as

vjn+1vjn=c9+ic10.(29)

The scheme will be stable if

|vjn+1vjn|2≤1.(30)

Which gives the stability conditions c92+c102≤1. If the scheme satisfies this condition, it will remain stable; otherwise, an unstable solution will be obtained. The stability condition is found for a scalar convection-diffusion problem.

To find the convergence condition for the system of convection diffusion equations, consider the following equation.

∂P∂t=A¯1∂P∂y+A¯2∂2P∂y2.(31)

where P is a vector and A¯1,A¯2 are matrices. The proposed scheme, when applied to Eq. (31), takes the form

P¯jn+1=12Pjne2Δt+12Pjne−Δt+(e2Δt+e−Δt−2){A¯1A1−1B1Pjn+A¯2A2−1B2Pjn−12Pjn},(32)

P=jn+1=aPjn+bP¯jn+1+14Δt{A¯1A1−1B1P¯jn+1+A¯2A2−1B2P¯jn+1},(33)

Pjn+1=cPjn+dP¯¯jn+1+23Δt{A¯1A1−1B1P¯¯jn+1+A¯2A2−1B2P¯¯jn+1}.(34)

Theorem: The proposed scheme (32)–(34) converges conditionally for the vectored matrix Eq. (31).

Proof: To prove this Theorem, consider the first stage of the exact scheme as:

P¯jn+1=12Pjne2Δt+12Pjne−Δt+(e2Δt+e−Δt−2){A¯1A1−1B1Pjn+A¯2A2−1B2Pjn−12Pjn}.(35)

Subtracting the first stage of the proposed and exact scheme, and let Pjn−pjn=Ejn, etc., it yields

E¯jn+1=12Ejne2Δt+12Ejne−Δt+(e2Δt+e−Δt−2){A¯1A1−1B1Ejn+A¯2A2−1B2Ejn−12Ejn}.(36)

Applying the norm ‖⋅‖∞ on both sides of Eq. (36) that yields

E¯n+1≤12Ene2Δt+12Ene−Δt+(e2Δt+e−Δt−2){‖A¯1A1−1B1‖∞En+‖A¯2A2−1B2‖∞En+12En}.(37)

Inequality (37) can be rewritten as

E¯n+1≤γ1En,(38)

where γ1=12e2Δt+12e−Δt+(e2Δt+e−Δt−2){‖A¯1A1−1B1‖∞+‖A¯2A2−1B2‖∞}.

The second stage of the exact scheme for Eq. (31) is written as

P=jn+1=aPjn+bP¯jn+1+14Δt{A¯1A1−1B1P¯jn+1+A¯2A2−1B2P¯jn+1}.(39)

The corresponding error equation is expressed as

E=jn+1=aEjn+bE¯jn+1+14Δt{A¯1A1−1B1E¯jn+1+A¯2A2−1B2E¯jn+1}.(40)

Applying the norm ‖⋅‖∞ on both sides of Eq. (40) as

E=n+1≤aEn+bE¯n+1+14Δt{‖A¯1A1−1B1‖∞E¯n+1+‖A¯2A2−1B2‖∞E¯n+1}.(41)

Rewrite inequality (41) as

E=n+1≤aEn+γ2E¯n+1.(42)

where γ2=b+14Δt{‖A¯1A1−1B1‖∞+‖A¯2A2−1B2‖∞}.

Substituting inequality (38) into inequality (42) results in

E=n+1≤γ3En,(43)

where γ3=a+γ1γ2.

The third stage of the exact scheme can be written as

Pjn+1=cPjn+dP¯¯jn+1+23Δt{A¯1A1−1B1P¯¯jn+1+A¯2A2−1B2P¯¯jn+1}.(44)

The error equation corresponding to Eq. (44) is

Ejn+1=cEjn+dE¯¯jn+1+23Δt{A¯1A1−1B1E¯¯jn+1+A¯2A2−1B2E¯¯jn+1}.(45)

Applying the norm ‖⋅‖∞ on both sides of Eq. (45) as

En+1≤cEn+dE=n+1+23Δt{‖A¯1A1−1B1‖∞+‖A¯2A2−1B2‖∞}E=n+1.(46)

Rewrite inequality (46) as

En+1≤cEn+γ4E=n+1+N(O((Δt)3,(Δy)6)),(47)

where γ4=d+23Δt{‖A¯1A1−1B1‖∞+‖A¯2A2−1B2‖∞}.

Using inequality (43) in inequality (47) as

En+1≤γ5En+N(O((Δt)3,(Δy)6)),(48)

where γ5=c+γ3γ4.

Let n=0 in inequality (48), which yields

E1≤γ5E0+N(O((Δt)3,(Δy)6)).(49)

Since E0=0, so inequality (49) becomes

E1≤N(O((Δt)3,(Δy)6)).(50)

Let n=1 in inequality (48), which results in

E2≤γ5E1+N(O((Δt)3,(Δy)6)).(51)

Using inequality (50) in inequality (51) gives

E2≤(1+γ5)+N(O((Δt)3,(Δy)6)).(52)

If this continued, then for finite n

En≤(1+γ5+…+γ5n−1)N(O((Δt)3,(Δy)6))=(1−γ5n1−γ)N(O((Δt)3,(Δy)6)).(53)

By taking limitas n→∞ the infinite series … is a finite geometric series that converges if |γ5|<1.□

4  Problem Statement

Consider one-dimensional, unsteady, laminar, incompressible, homogeneous in the x∗-direction (i.e., ∂∂x=0) Eyring-Prandtl nano fluid flow over Riga plate. The Riga plate comprises an array of electrodes and magnets that generate a steady Lorentz force acting tangentially to the fluid. Let x∗-axis is placed vertically and y∗-axis is perpendicular to it. The Riga plate lies along x∗-axis parallel to x∗-axis. The fluid occupies the region y∗>0. Let the plate suddenly move along the positive x∗-axis. A slip boundary condition is imposed at the wall. Let Tw and Cw are temperature and nanoparticle concentration at the plate. Ambient temperature and nanoparticle concentration are far away from the plate and are denoted, respectively, by T∞ and C∞.

Model Assumptions and Coordinate System: A Cartesian coordinate system is used: x∗-axis is aligned vertically along the plate (direction of flow). y∗-axis is normal to the plate, pointing into the fluid domain. The flow is assumed to be: Unsteady, i.e., velocity, temperature, and concentration vary with time. Incompressible, meaning the fluid density is constant. Laminar and one-dimensional in the y∗-direction; all ∂/∂x=0. Homogeneous and Newtonian in bulk, but governed by the Eyring-Prandtl non-Newtonian viscosity model. The fluid is: A nanofluid, i.e., a base fluid containing suspended nanoparticles. Electrically conducting and influenced by electromagnetic forces generated by the Riga plate. The Riga plate is: A surface embedded with alternating magnets and electrodes that produce a steady Lorentz force along the x∗-direction. Subject to a velocity slip boundary condition, rather than a no-slip. Heat and mass transfer considerations include: Thermophoresis and Brownian motion effects (key features of nanofluid dynamics). Viscous dissipation due to Eyring–Prandtl fluid behavior. A chemical reaction term acting as a sink in the energy and concentration equations. The ambient conditions are: Far from the plate, fluid is at rest, and at ambient temperature T∞ and concentration C∞. Additional assumptions: No external pressure gradient is imposed. Constant thermal conductivity and mass diffusivity. The governing equations are written under dimensional and dimensionless forms to facilitate analysis. The governing equations are expressed as:

∂u∗∂t∗=β¯C¯(∂2u∗∂y∗21+(1C¯∂u∗∂y∗)2)+gβT(T−T∞)+gβC(C−C∞)+πj∘M∘8ρexp⁡(−πy∗a1),(54)

∂T∂t∗=α∂2T∂y∗2+β¯ρCpsinh−1⁡(1C¯∂u∗∂y∗)∂u∗∂y∗+τ(DB∂C∂y∗∂T∂y∗+DTT∞(∂T∂y∗)2)−k1(C−C∞),(55)

∂C∂t∗=DB∂2C∂y∗2+DTT∞∂2T∂y∗2.(56)

Momentum Eq. (54) (x∗-direction): ∂u∗∂t∗ represents an unsteady change in velocity; it tracks how the fluid accelerates or decelerates over time. ∂2u∗∂y∗2 is the Momentum diffusion (viscous force) in y∗-direction, it captures the spread of velocity due to internal friction. 1+(1C¯∂u∗∂y∗)2 represents the Eyring–Prandtl nonlinearity (shear-thinning), it models how viscosity varies with shear, common in polymers, biofluids, etc. β¯C¯ is the stress scale from the Eyring–Prandtl model controls how strongly the fluid resists deformation. gβT(T−T∞) and gβC(C−C∞) are the thermal and Solutal buoyancy forces, respectively. πj∘M∘8ρe−(πa1)y∗ is the Lorentz force (electromagnetic) from the Riga plate, it drives the fluid forward along the wall with exponentially decaying intensity. The Eyring–Prandtl term captures non-Newtonian shear-thinning behavior in smart fluids, slurries, and synthetic lubricants. Buoyancy terms are important in free-convection nanofluid flows, e.g., in cooling devices or bioconvection, and the Riga plate force is the main driver of the flow, representing active electromagnetic pumping used in microfluidics, MHD accelerators, or coating technologies.

Energy Eq. (55): ∂T∂t∗ represents the rate of change of temperature, also referred to as unsteady heat conduction, α∂2T∂y∗2 is the Thermal diffusion spreads heat within the fluid, β¯ρCpsinh−1⁡(1C¯∂u∗∂y∗)∂u∗∂y∗ is the viscous dissipation (frictional heating), it converts mechanical work into thermal energy, τDB∂C∂y∗∂T∂y∗ is the Brownian-motion-induced coupling, which affects heat due to nanoparticle concentration gradients, τDTT∞(∂T∂y∗)2 is the Thermophoresis effect, when nanoparticles migrate from hot to cold regions, affecting the thermal field and −k1(C−C∞) represents the Energy sink due to a chemical reaction; it removes heat from exothermic or endothermic processes.

Nanoparticle Concentration Eq. (56): ∂C∂t∗ represents the rate of change of nanoparticle concentration, it captures nanoparticle transport over time, DB∂2C∂y∗2 is the Mass diffusion (Brownian motion), in a random motion of particles, increases the spread and DTT∞∂2T∂y∗2 Thermophoresis (temperature-driven diffusion) causes particles to move from hot to cold areas.

Subject to the following initial and boundary conditions

u∗=0,T=0,C=0whent∗=0u∗=uw+Ls∂u∗∂y∗,T=Tw,C=Cwaty∗=0u∗→0,T→T∞,C→C∞wheny∗→∞}.(57)

Physical Significance of initial and Boundary Conditions (Dimensional Form): Initial condition at time t∗=0, the fluid is at rest. It is at ambient temperature and ambient nanoparticle concentration. Boundary conditions at the wall y∗=0, unlike the no-slip condition in classical flows, here: u∗(t∗,0)=uw+Ls∂u∗∂y∗|y∗=0, uw is the reference wall velocity (often a constant or specified function) and Ls is the slip length, indicating how much fluid slides along the wall. T=Tw,C=Cw is the Dirichlet condition prescribes heat and nanoparticle injection. Far-field conditions as y∗→∞, far from the plate: The fluid is undisturbed (no motion). It maintains ambient temperature and ambient nanoparticle concentration.

Where u∗ is the horizontal component of velocity, C¯ is the rate parameter, β¯ is a stress parameter, M∘ is a magnetic field strength, j∘ is current density, a1 width of electrodes, g is gravity, βT is thermal Expansion, βC is solutal Expansion, α is thermal diffusivity, DB is the Brownian motion coefficient, DT represents the thermophoresis coefficients, cp is the specific heat capacity, k1 is the reaction rate parameter.