Dynamic simulation of non-Newtonian boundary layer flow: An enhanced exponential time integrator approach with spatially and temporally variable heat sources

Muhammad Shoaib Arif; Kamaleldin Abodayeh; Yasir Nawaz

doi:10.1515/phys-2024-0034

Artikel Open Access

Dynamic simulation of non-Newtonian boundary layer flow: An enhanced exponential time integrator approach with spatially and temporally variable heat sources

Muhammad Shoaib Arif , Kamaleldin Abodayeh und Yasir Nawaz

Veröffentlicht/Copyright: 30. Mai 2024

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Open Physics Band 22 Heft 1

Abstract

Scientific inquiry into effective numerical methods for modelling complex physical processes has led to the investigation of fluid dynamics, mainly when non-Newtonian properties and complex heat sources are involved. This paper presents an enhanced exponential time integrator approach to dynamically simulate non-Newtonian boundary layer flow with spatially and temporally varying heat sources. We propose an explicit scheme with second-order accuracy in time, demonstrated to be stable through Fourier series analysis, for solving time-dependent partial differential equations (PDEs). Utilizing this scheme, we construct and solve dimensionless PDEs representing the flow of Williamson fluid under the influence of space- and temperature-dependent heat sources. The scheme discretizes the continuity equation of incompressible fluid and Navier–Stokes, energy, and concentration equations using the central difference in space. Our analysis illuminates how factors affect velocity, temperature, and concentration profiles. Specifically, we observe a rise in temperature profile with enhanced coefficients of space and temperature terms in the heat source. Non-Newtonian behaviours and geographical/temporal variations in heat sources are critical factors influencing overall dynamics. The novelty of our work lies in developing an explicit exponential integrator approach, offering stability and second-order accuracy, for solving time-dependent PDEs in non-Newtonian boundary layer flow with variable heat sources. Our results provide valuable quantitative insights for understanding and controlling complex fluid dynamics phenomena. By addressing these challenges, our study advances numerical techniques for modelling real-world systems with implications for various engineering and scientific applications.

Keywords: exponential integrator; stability; convergence; Williamson fluid; heat source

1 Introduction

A great application of boundary layer flows with non-Newtonian liquids is important to study as it has valuable practical and theoretical applications in various natural and industrial environments. The complications in the flow dynamics of non-Newtonian fluids distinct from that of Newtonian fluids are mainly because of their complicated rheological features. Williamson liquid also belongs to the class of non-Newtonian fluids. To gain an insight into the complex nature of the Williamson liquid boundary layer flow process, a modified exponential integrator method was used for this project. The flow characteristics are also considerably different when a heat source that varies with both space and temperature is taken into account, adding another layer of subtlety for a more practical analysis for the study of the boundary layer flow of such a fluid.

Numerous biological and commercial processes use a Williamson fluid model, which is an example of a non-Newtonian fluid. The Williamson fluids have viscoelastic properties and also have some unique rheological characteristics, such as shear-thinning or shear-thickening behaviour. These unusual properties of Williamson fluids require them to be modelled using a proper mathematical model. However, traditional numerical methods do not ensure an accurate simulation of the non-Newtonian flows. In the investigation of non-Newtonian flows, our research utilizes a new technique, that is, a modified exponential integrator approach to eliminate the disadvantages of the traditional numerical simulation method. Solving stiff ordinary differential equations with non-Newtonian fluid models has its separate difficulties. As a result, this study employs a particular numerical method to eliminate the problem.

The restricted application of Newtonian fluids has sparked interest in studying non-Newtonian fluids. Some non-Newtonian fluids are starch, lubricating sprays, honey, ketchup, and countless more. The rheological characteristics of fluids were demonstrated by Williamson in 1929 using a non-Newtonian model [1]. It is known as the Williamson fluid model in academic publications. Nadeem and Hussain [2] studied the effect of heat transfer through a magnetized Williamson fluid across a continuously stretched sheet. Using the Cattaneo-Christov (CC) heat flow model, Amjad et al. [3] investigated the magnetohydrodynamic (MHD) Williamson nanofluid moving across an exponential stretching sheet. The peristaltic transport of a Williamson nanofluid in a tapered asymmetric channel under the effect of a thermal radiation parameter is theoretically investigated by Kothandapani and Prakash [4]. Taking into account two scenarios for heat transfer, the prescribed exponential order surface temperature and the prescribed exponential order heat flux, Ahmed and Akbar [5] detailed the analysis of an MHD flow of Williamson nanofluid across an exponentially porous stretched surface. Using models and analyses, Hayat et al. [6] sought to determine how well magnetic nanoparticles and melting heat transfer performed in the stretched flow. Abbas et al. [7] studied the computational assessment of the impact of MHDs and changing density on Williamson Sakiadis flow in a porous medium. Shafiq et al. [8] presented a new approach to solving numerical problems using artificial neural networks that combines the Levenberg–Marquard algorithm with multi-layer perceptron feed-forward back-propagation. This approach can understand radiation, heat generation/absorption, and unsteady electrically conducting Williamson liquid flow along porous stretching surfaces. Ahmed et al. [9] examined the mechanics of mixed-heterogeneous two-dimensional Williamson fluid flows across a nonlinear extending curved surface while considering convective boundary conditions and homogeneous–heterogeneous responses. Ahmed et al. [10] focused on analysing the heat flux mechanism in a MHD mixed convective flow of a Williamson-type fluid over an exponential stretching porous curved surface.

The esteemed Nobel laureate Hannes Alfven is believed to have proposed the theory of MHD. It explores the relationship between fluid mechanics and electromagnetism, illustrating the impact of a magnetic field on a conductive fluid. Some of the applications of MHDs involve magnetic endoscopy, tissue temperature analysis, centrifugal pumps, cell separation, blood flow, and cancer tumour therapy. An investigation was conducted on nanofluids’ two-dimensional boundary layer flow undergoing a chemical reaction with thermal radiation. The medium of study was a nonlinear stretched plate [11]. Turkyilmazoglu [12] utilized an analytical technique to handle the MHD flow. A study by Mabood et al. [13] investigated the flow of magnetic fluid along a vertical axis while the fluid was spinning. Talihnoee et al. [14] studied nanofluids to examine their potential for generating MHD entropy and spontaneous convection. In the study conducted by Rashidi et al. [15], it was proposed that the stretching of porous sheets could be influenced by heat radiation and thermodiffusion in the MHD Williamson fluid. Madhu et al. [16] studied the minimization of entropy generation through the thermodynamic second law, where we consider heat transfer and non-Newtonian fluid (Williamson fluid) flow via a micro-channel. Mishra et al. [17] studied Williamson fluid’s micro-rotational micro-hydrodynamic flow through a non-Darcy porous material. Almaneea [18] investigated the impact of hybrid nanoparticles on heat and mass transport processes in both homogeneous and heterogeneous chemical reactions. The finite-element method (FEM) is utilized to solve models numerically. Reddy et al. [19] investigated the mass and heat transfer properties of MHD Williamson nanofluids in porous media with CC double diffusion on a stretching surface. Asjad et al. [20] examined the impact of Brownian motion and thermophoresis diffusion in the flow of a non-Newtonian Williamson fluid over an exponentially stretched sheet, considering the impacts of thermal radiation and the bioconvection of microorganisms.

A higher stretching ratio reduces the friction factor. It increases the heat transfer rate, according to research by Raju et al. [21], on the impact of a space and temperature-dependent heat source on the flow of a nanofluid in a boundary layer across a uniformly thick surface when a nonlinearly permeable stretching sheet is used. Mather initially documented a viscous liquid’s flow across a paraboloid surface in 1961 [22]. After that, Lee [23] discussed the flow of an incompressible fluid over a thin needle, precisely a non-uniformly thick item, and provided the boundary layer equation that governs this motion. Applying the idea to boundary layer flow on thin paraboloids, Cebeci et al. [24] expanded upon it. Upstream solutions for continuous, viscous flow via a paraboloid and the boundary layer on a paraboloid of revolution were studied by Miller [25,26]. Veldman [27] proposed a numerical method for solving the Navier–Stokes equations in the context of flow around a paraboloid of revolution.

In fluid dynamics, fluids can be classified as Newtonian or non-Newtonian. The viscous strains caused by Newtonian fluid flow are clearly proportionate to the local strain rate. For non-Newtonian fluids, the viscous stresses caused by the flow do not relate linearly and could even change with time. Three main types of fluids are not Newtonian: viscoelastic, time-dependent, and time-independent. Temperature, shear rate, and time all play a role in determining the viscosity of non-Newtonian fluid time-dependent. Paint and yogurt are examples of thixotropic materials that decrease viscosity with time. At the same time, gypsum paste is a rheotactic material that increases in viscosity with time. In contrast, the viscosity of time-independent non-Newtonian fluids is solely affected by temperature and shear rate.

Shear thickening, in contrast, is a dilatant process where the viscosity increases in proportion to the shear rate. Lastly, plasticity is the property that states flow can only occur under very particular shear stress conditions. The shear rate at a specific place in a time-independent non-Newtonian fluid is solely determined by the stress at that point.

A model equation for describing the flow of pseudoplastic fluids was proposed and experimentally validated by Williamson [1], who also addressed the flow of pseudoplastic materials. Williamson fluid is an example of visco-inelastic fluids from a rheological perspective. According to Nadeem et al. [28], Williamson fluid is a non-Newtonian fluid with a shear thinning property, meaning that its viscosity decreases as the rate of shear stress increases. Specialists have documented fluid flows on different surfaces since appropriate models for Casson and Williamson were developed.

Heat is the amount of energy transferred between two bodies with different magnitudes of temperature, whereas temperature is the degree to which an object or substance is hot or cold. Thermal expansion is described as the propensity of matter to change its shape, area, and volume in response to a substantial alteration in temperature caused by heat transfer.

Another way to look at heat is as an energy source that moves freely between an object and its environment. Three primary pathways exist for heat transfer: conduction, radiation, and convection. The term “convection” describes the movement of heat across a fluid medium when temperature gradients cause density variations in the fluid. A process known as free or natural convection takes place rather than the presence of an external heat source. Fluids near a heat source experience free convection, which causes them to lose density and ascend. The surrounding colder fluid subsequently replaces it.

Shehzad et al. [29] examined the scenario of heat and mass convection at the surface after considering this theory and considering flow induced by a bidirectional stretched surface. Heat transmission in an unstable, convection nanofluid with gyrotactic microbes and nanoparticles was documented by Motsa and Animasaun [30] using the paired quasi-linearization method. Abbasi et al. [31] investigated Maxwell nanofluid mixed convection flow. The first unstable stage, when the temperature increases the maximum, is when the convective acceleration element in the energy equation has a minimal impact, according to Sandeep et al. [32], who analysed three-dimensional Casson fluid flow over a zero-temperature surface. Bhatti et al. [33] provided a detailed discussion of the presence of diamond ( C ) and silica ( Si O 2 ) nanoparticles in a water-based hybrid nanofluid floating on an exponentially elastic surface. In Bhatti et al. [34], a mathematical model for coating boundary layer flow of a Maxwell viscoelastic fluid with non-Fourier heat flux and nonlinear quadratic convection is proposed. Heat source/sink effects and nonlinear quadratic thermal radiation are also considered. Lie symmetry transformations are used. Nadeem et al. [35] analysed entropy creation by investigating irreversible sources in the steady flow of a non-Newtonian Williamson fluid.

1.1 Novelty of this study

The study mentioned above stands out in numerical fluid dynamics and heat transfer methods due to its numerous new contributions, particularly in the context of non-Newtonian boundary layer flow affected by spatial and temperature-dependent heat sources. In a nutshell, the following are the main innovations of this study:

Exponential integrator scheme: The main novelty is a new method for solving time-dependent partial differential equations (PDE) using an exponential integrator technique. This scheme provides a novel way to handle the problems given by dynamic systems effectively; it is explicit and has second-order temporal precision.
Stability analysis: Using Fourier series analysis, the research thoroughly examines the stability of the suggested exponential integrator system. To guarantee the correctness and dependability of the results, it is essential to establish the stability of a numerical scheme; the offered analysis adds rigor to the investigation.
Application to non-Newtonian fluid dynamics: The suggested exponential integrator is put to new use in modelling Williamson fluid flow, an innovative application of the numerical technique. Understanding how non-Newtonian fluids differ from conventional Newtonian fluid difficult endeavour, but studying their dynamics still adds complexity.
Space- and temperature-dependent heat source: This research focuses on analysing the effects of a heat source with varying coefficients in both space and time. This characteristic enhances the mathematical model’s realism in scenarios involving dynamic and diverse heat sources. The study stands out because it investigates how these changes affect the system’s behaviour.
Comprehensive discretization approach: To provide a thorough numerical framework, the research utilized a discretization strategy that included a first-order scheme for the continuity equation as well as the proposed exponential integrator for the Navier–Stokes equation, energy equation, and concentration equation with central difference in space. This comprehensive strategy enhances the effectiveness and precision of the numerical solution.

Our research aims to address the following primary question: How can we design efficient numerical models to capture the complex structure of non-Newtonian boundary layer flow when heat sources are both geographically and temporally variable? Reliable numerical techniques for solving time-dependent PDEs regulating such processes are the focus of our work, guided by this overall concern.

Our investigation goes beyond the standard boundaries by including space- and temperature-dependent heat sources in the analysis. By including these extra factors, we hope to replicate more realistic situations that users may face in real-world applications. Because of their centrality in many environmental and engineering situations, heat sources significantly influence the distribution of temperatures and, by extension, the dynamics of fluids. By including space- and temperature-dependent heat sources, the study is improved, the underlying physical processes are better represented, and the study’s applicability is widened.

First, a complete mathematical framework for modelling non-Newtonian Williamson fluid boundary layer flow using the modified exponential integrator for improved numerical accuracy is introduced; second, the effect of space- and temperature-dependent heat sources on the flow properties is investigated; and finally, conclusions are drawn from these results. Our goal in this systematic investigation is to illuminate the complex relationship between non-Newtonian fluid dynamics, numerical methods, rheological parameters, and heat sources. Understanding complicated fluid behaviours is crucial for engineers, environmental scientists, and others to optimize their processes.

Mathematical boundary layer flow models are expressed as PDEs. Further, these PDEs are reduced into dimensionless partial or ordinary differential equations. Matlab solver, analytical methods, and finite or FEMs can solve ordinary differential equations. In some situations, the accuracy of the solution is only required instead of any specific method. Among the mentioned methods, the finite difference method can be used to solve these equations. In this contribution, an exponential integrator is proposed to solve the dimensionless set of PDEs arising in the flow phenomenon of flow over the flat and oscillatory plates. The method is clearly defined and achieves second-order accuracy in time. Spatial discretization is performed using second-order central formulas.

2 Exponential time integrator

A two-stage explicit scheme is proposed that can be used to solve time-dependent problems in science and engineering. The scheme discretizes time-dependent terms in a given PDE. Any other scheme can be used for space discretization. To propose a time integrator, consider an equation of the form [36]

(1) ∂ v ∂ t = G v , ∂ v ∂ x , ∂ v ∂ y , ∂ 2 v ∂ y 2 .

Let Eq. (1) be rearranged as

(2) ∂ v ∂ t = v + F ,

where F = G − v .

The initial step of the scheme is described as

(3) v ̅ i , j n + 1 = v i , j n e Δ t + ( e Δ t − 1 ) ∂ v ∂ t i , j n − v i , j n ,

where Δ t is the time step size; this stage predicts the solution, and the solution will be corrected in the corrector stage. The corrector stage of the scheme is expressed as follows:

(4) v i , j n + 1 = 1 2 ( v i , j n e Δ t + v ̅ i , j n + 1 ) + ( e Δ t − 1 ) { a ( G i , j n − c v i , j n ) + b ( G ̅ i , j n + 1 − v ̅ i , j n + 1 ) } ,

where G ̅ i , j n + 1 = G v ̅ i , j n + 1 , ∂ v ̅ ∂ x i , j n + 1 , ∂ v ̅ ∂ y i , j n + 1 ∂ 2 v ̅ ∂ y 2 i , j n + 1 .

The second stage (4) contains three parameters a , b , and c . The values of these parameters will be determined by utilizing the Taylor series expansion. To do this, re-write Eq. (4) as

(5) v i , j n + 1 = 1 2 ( v i , j n e Δ t + v ̅ i , j n + 1 ) + ( e Δ t − 1 ) a ∂ v ∂ t i , j n − c v i , j n + b ∂ v ̅ ∂ t i , j n + 1 − v ̅ i , j n + 1 .

The Taylor series expansion for v i , j n + 1 is given as

(6) v i , j n + 1 = v i , j n + Δ t ∂ v ∂ t i , j n + ( Δ t ) 2 ∂ 2 v ∂ t 2 i , j n + O ( ( Δ t ) 3 ) .

Substituting the Taylor series expansion (6) into Eq. (5), the following is obtained:

(7) v i , j n + Δ t ∂ v ∂ t i , j n + ( Δ t ) 2 ∂ 2 v ∂ t 2 i , j n = 1 2 ( v i , j n e Δ t + v ̅ i , j n + 1 ) + ( e Δ t − 1 ) a ∂ v ∂ t i , j n − c v i , j n + b ∂ v ̅ ∂ t i , j n + 1 − v ̅ i , j n + 1 .

When the first stage of the scheme (3) is substituted into Eq. (7), the result is

(8) v i , j n + Δ t ∂ v ∂ t i , j n + ( Δ t ) 2 ∂ 2 v ∂ t 2 i , j n = v i , j n e Δ t + e Δ t − 1 2 ∂ v ∂ t i , j n − v i , j n + ( e Δ t − 1 ) a ∂ v ∂ t i , j n − c v i , j n + b e Δ t ∂ v ∂ t i , j n + ( e Δ t − 1 ) ∂ 2 v ∂ t 2 i , j n − ∂ v ∂ t i , j n − v i , j n e Δ t − ( e Δ t − 1 ) ∂ v ∂ t i , j n − v i , j n .

Equating the coefficients of v i , j n , ∂ v ∂ t i , j n and ∂ 2 v ∂ t 2 i , j n on both sides of Eq. (8) yields

(9) 1 = e Δ t − e Δ t − 1 2 − ac ( e Δ t − 1 ) − b e Δ t ( e Δ t − 1 ) + b ( e Δ t − 1 ) 2 ,

(10) Δ t = e Δ t − 1 2 + a ( e Δ t − 1 ) + b e Δ t ( e Δ t − 1 ) − b ( e Δ t − 1 ) 2 − b ( e Δ t − 1 ) 2 ,

(11) ( Δ t ) 2 2 = b ( e Δ t − 1 ) 2 .

Upon solving Eqs. (9)–(11), the values of parameters a , b , and c are given as follows:

(12) a = 1 + 2 Δ t + 2 ( Δ t ) 2 − 2 e Δ t − 2 Δ t e Δ t − ( Δ t ) 2 e Δ t + e 2 Δ t 2 ( e Δ t − 1 ) 2 b = ( Δ t ) 2 2 ( e Δ t − 1 ) 2 c = 1 − ( Δ t ) 2 − 2 e Δ t + e 2 Δ t 1 + 2 ( Δ t ) 2 − 2 e Δ t − ( Δ t ) 2 e Δ t + e 2 Δ t + 2 Δ t − 2 Δ t e Δ t .

Let G = β 1 ∂ v ∂ x + β 2 ∂ v ∂ y + β ∂ 2 v ∂ y 2 , then the time discretization for Eq. (1) is

(13) v ̅ i , j n + 1 = v i , j n e Δ t + ( e Δ t − 1 ) β 1 ∂ v ∂ x i , j n + β 2 ∂ v ∂ y i , j n + β 3 ∂ 2 v ∂ y 2 i , j n − v i , j n ,

(14) v i , j n + 1 = 1 2 ( v i , j n e Δ t + v ̅ i , j n + 1 ) + ( e Δ t − 1 ) a β 1 ∂ v ∂ x i , j n + β 2 ∂ v ∂ y i , j n + β 3 ∂ 2 v ∂ y 2 i , j n − c v i , j n + b β 1 ∂ v ̅ ∂ x i , j n + 1 + β 2 ∂ v ̅ ∂ y i , j n + 1 + β 3 ∂ 2 v ̅ ∂ y 2 i , j n + 1 − v ̅ i , j n + 1 .

Consider the second-order central difference formula for space discretization, then Eqs. (13) and (14) can be expressed as follows:

(15) v ̅ i , j n + 1 = v i , j n e Δ t + ( e Δ t − 1 ) { β 1 δ x v i , j n + β 2 δ y v i , j n + β 3 δ y 2 v i , j n − v i , j n } ,

(16) v i , j n + 1 = 1 2 ( v i , j n e Δ t + v ̅ i , j n + 1 ) + ( e Δ t − 1 ) { a ( β 1 δ x v i , j n + β 2 δ y v i , j n + β 3 δ y 2 v i , j n − c v i , j n ) + b ( β 1 δ x v ̅ i , j n + 1 + β 2 δ y v ̅ i , j n + 1 + β 3 δ y 2 v ̅ i , j n + 1 − v ̅ i , j n + 1 ) } ,

where δ x v i , j n = v i + 1 , j n − v i − 1 , j n 2 ( Δ x ) , δ y v i , j n = v i , j + 1 n − v i , j − 1 n 2 ( Δ y ) , δ y 2 v i , j n = v i , j + 1 n − 2 v i , j n + v i , j − 1 n ( Δ y ) 2 and Δ x , Δ y are spatial step sizes in the x and y directions, respectively.

3 Stability analysis

The stability analysis of finite difference methods for linear differential equations is found in the literature by applying Von Neumann stability analysis. It can also be used for finding stability conditions of nonlinear differential equations discretized by finite difference methods. Equations are linearized to apply this analysis for finite difference methods on nonlinear differential equations. So, in this, just an estimation of differential equations will be given. To apply this method, the dependent variable in the difference equation is transformed to Fourier series. Instead of considering the whole series, only one component is employed for transformation. Consider the following transformation to initiate the application procedure for these criteria:

(17) v ̅ i , j n + 1 = P ̅ n + 1 e iI ψ 1 e jI ψ 2 , v i , j n = P n e iI ψ 1 e jI ψ 2 v ̅ i ± , j n + 1 = P ̅ n + 1 e ( i ± 1 ) I ψ 1 e jI ψ 2 , v ̅ i , j ± 1 n + 1 = P ̅ n + 1 e iI ψ 1 e ( j ± 1 ) I ψ 2 v i ± 1 , j n = P n e ( i ± 1 ) iI ψ 1 e jI ψ 2 , v i , j ± 1 n = P n e iI ψ 1 e ( j ± 1 ) I ψ 2 v i , j n + 1 = P n + 1 e iI ψ 1 e jI ψ 2 .

Utilizing some transformations from (17) into the first stage of the scheme (15) yields

(18) P ̅ n + 1 e iI ψ 1 e jI ψ 2 = P n e iI ψ 1 e jI ψ 2 e Δ t + ( e Δ t − 1 ) × β 1 ( e ( i + 1 ) I ψ 1 e jI ψ 2 − e ( i − 1 ) I ψ 1 e jI ψ 2 ) 2 ( Δ x ) P n + β 2 ( e iI ψ 1 e ( j + 1 ) I ψ 2 − e iI ψ 1 e ( j − 1 ) I ψ 2 ) 2 ( Δ y ) P n + β 3 ( e iiI ψ 1 e ( j + 1 ) I ψ 2 − 2 e iI ψ 1 e jI ψ 2 + e iI ψ 1 e ( j − 1 ) I ψ 2 ) ( Δ y ) 2 P n − P n e iI ψ 1 e jI ψ 2 .

Now, dividing both sides of Eq. (18) by e iI ψ 1 e jI ψ 2 gives

(19) P ̅ n + 1 = P n e Δ t + ( e Δ t − 1 ) β 1 ( e I ψ 1 − e − I ψ 1 ) 2 ( Δ x ) P n + β 2 ( e I ψ 2 − e − I ψ 2 ) 2 ( Δ y ) P n + β 3 ( e − I ψ 2 − 2 + e − I ψ 2 ) ( Δ y ) 2 P n − P n .

By using De Movier’s theorem, Eq. (19) is expressed as follows:

(20) P ̅ n + 1 = P n e Δ t + ( e Δ t − 1 ) β 1 2 I sin ψ 1 2 ( Δ x ) P n + β 2 2 I sin ψ 2 2 ( Δ y ) P n + β 3 ( 2 cos ψ 2 − 2 ) ( Δ y ) 2 P n − P n .

Re-write Eq. (20) as follows:

(21) P ̅ n + 1 = [ e Δ t + { c 1 I sin ψ 1 + c 2 I sin ψ 2 + 2 d 1 ( cos ψ 2 − 1 ) − ( e Δ t − 1 ) } ] P n ,

where

c 1 = β 1 ( e Δ t − 1 ) Δ x , c 2 = β 2 ( e Δ t − 1 ) Δ y , d 1 = β 3 ( e Δ t − 1 ) ( Δ y ) 2 .

Similarly, utilizing transformations from Eq. (17) into the second stage of the scheme (16) yields

(22) P n + 1 = 1 2 ( P n e Δ t + P ̅ n + 1 ) + ( e Δ t − 1 ) a β 1 I sin ψ 1 ( Δ x ) P n + β 2 I sin ψ 2 ( Δ y ) P n + 2 β 3 ( cos ψ 2 − 1 ) ( Δ y ) 2 P n − c P n + b β 1 I sin ψ 1 ( Δ x ) P ̅ n + 1 + β 2 I sin ψ 2 ( Δ y ) P ̅ n + 1 + 2 β 3 ( cos ψ 2 − 1 ) ( Δ y ) 2 P ̅ n + 1 − P ̅ n + 1 .

Eq. (22) is re-written as follows:

(23) P n + 1 = 1 2 e Δ t + a c 1 I sin ψ 1 + a c 2 sin ψ 2 + 2 a d 1 ( cos ψ 2 − 1 ) − c ( e Δ t − 1 ) P n + 1 2 + b c 1 I sin ψ 1 + b c 2 I sin ψ 2 + 2 b d 1 ( cos ψ 2 − 1 ) − ( e Δ t − 1 ) P ̅ n + 1 .

By using Eq. (21) in Eq. (23), we obtain

(24) P n + 1 = 1 2 e Δ t + a c 1 I sin ψ 1 + a c 2 sin ψ 2 + 2 a d 1 ( cos ψ 2 − 1 ) − c ( e Δ t − 1 ) P n + 1 2 + b c 1 I sin ψ 1 + b c 2 I sin ψ 2 + 2 b d 1 ( cos ψ 2 − 1 ) − ( e Δ t − 1 ) ( e Δ t + c 1 I sin ψ 1 + c 2 I sin ψ 2 + 2 d 1 ( cos ψ 2 − 1 ) − ( e Δ t − 1 ) ) P n .

The stability condition is given as

(25) P n + 1 P n 2 = 1 2 e Δ t − c ( e Δ t − 1 ) + a Re ( z ) + 1 2 − ( e Δ t − 1 ) + b Re ( z ) Re ( z ) − b Im ( z ) 2 + ( a Im ( z ) + Im ( z ) ) 2 < 1 ,

where

z = c 1 I sin ψ 1 + c 2 sin ψ 2 + 2 a d 1 ( cos ψ 2 − 1 )

Inequality (25) gives the stability condition of the proposed finite difference method for the linear scalar convection–diffusion equation. Further, this contribution convergence analysis for the system of convection–diffusion equation is also provided. For this purpose, consider the system of convection–diffusion equation

(26) ∂ u ∂ t = A ∂ u ∂ x + B ∂ u ∂ y + C 1 ∂ 2 u ∂ y 2 + Du .

By employing the proposed scheme in Eq. (26), the discretized equation can be expressed as

(27) u ̅ i , j n + 1 = u i , j n e Δ t + ( e Δ t − 1 ) { A δ x u i , j n + B δ y u i , j n + C 1 δ y 2 u i , j n + D u i , j n } ,

(28) u i , j n + 1 = 1 2 ( u i , j n e Δ t + u ̅ i , j n + 1 ) + ( e Δ t − 1 ) { a ( A δ x u i , j n + B δ y u i , j n + C 1 δ y 2 u i , j n + D u i , j n − c u i , j n ) + b ( A δ x u ̅ i , j n + 1 + B δ y u ̅ i , j n + 1 + C 1 δ y 2 u ̅ i , j n + 1 + D u ̅ i , j n + 1 − u ̅ i , j n + 1 ) } .

Theorem 1

The proposed exponential integrator schemes (27) and (28) converge for the system of Eq. (26).

Proof

For proving this theorem, consider the exact scheme for Eq. (26)□

(29) U ̅ i , j n + 1 = U i , j n e Δ t + ( e Δ t − 1 ) { A δ x U i , j n + B δ y U i , j n + C 1 δ y 2 U i , j n + D U i , j n } ,

(30) U i , j n + 1 = 1 2 ( U i , j n e Δ t + U ̅ i , j n + 1 ) + ( e Δ t − 1 ) { a ( A δ x U i , j n + B δ y U i , j n + C 1 δ y 2 U i , j n + D U i , j n − c U i , j n ) + b ( A δ x U ̅ i , j n + 1 + B δ y U ̅ i , j n + 1 + C 1 δ y 2 U ̅ i , j n + 1 + D U ̅ i , j n + 1 − U ̅ i , j n + 1 ) } .

By subtracting the first stage of the proposed scheme (27) from the first stage of the exact scheme (29) and also let

U ̅ i , j n + 1 − u ̅ i , j n + 1 = e ̅ i , j n + 1 , U i , j n − u i , j n = e i , j n , U i ± 1 , j n − u i ± 1 , j n = e i ± 1 , j n , U i , j ± 1 n − u i , j ± 1 n = e i , j ± 1 n .

Then, the error equation is written as

(31) e ̅ i , j n + 1 = e i , j n e Δ t + ( e Δ t − 1 ) { A δ x e i , j n + B δ y e i , j n + C 1 δ y 2 e i , j n + D e i , j n } .

Apply norm ‖ . ‖ on both sides of Eq. (31) as follows:

(32) e ̅ 1 n + 1 ≤ e 1 n e Δ t + ( e Δ t − 1 ) ‖ A ‖ ∞ e 1 n Δ x + ‖ B ‖ ∞ e 1 n Δ y + ‖ C 1 ‖ ∞ 4 e 1 n ( Δ y ) 2 + ‖ D ‖ ∞ e 1 n .

Let

c ̅ 1 = ‖ A ‖ ∞ ( e Δ t − 1 ) Δ x , c ̅ 2 = ‖ B ‖ ∞ ( e Δ t − 1 ) Δ y , d ̅ 1 = ‖ C 1 ‖ ∞ ( e Δ t − 1 ) ( Δ y ) 2

Re-write inequality (32) as

(33) e ̅ 1 n + 1 ≤ e 1 n e Δ t + c ̅ 1 e 1 n + c ̅ 2 e 1 n + 4 d ̅ 1 e 1 n + ‖ D ‖ ∞ ( e Δ t − 1 ) e 1 n } .

Similarly, subtracting the second stage of the proposed scheme (28) from the first stage of the exact scheme (30) and let U i , j n + 1 − u i , j n + 1 = e i , j n + 1 , U ̅ i ± 1 , j n + 1 − u ̅ i ± 1 , j n + 1 = e ̅ i ± 1 , j n + 1 gives

(34) e i , j n + 1 = 1 2 ( e i , j n e Δ t + e ̅ i , j n + 1 ) + ( e Δ t − 1 ) { a ( A δ x e i , j n + B δ y e i , j n + C 1 δ y 2 e i , j n + D e i , j n − c e i , j n ) + b ( A δ x e ̅ i , j n + 1 + B δ y e ̅ i , j n + 1 + C 1 δ y 2 e ̅ i , j n + 1 + D e ̅ i , j n + 1 − e ̅ i , j n + 1 ) } .

Applying norm ‖ . ‖ on both sides of Eq. (34) gives

(35) e 1 n + 1 ≤ 1 2 ( e 1 n e Δ t + e ̅ 1 n + 1 ) + ( e Δ t − 1 ) | a | ‖ A ‖ ∞ e 1 n Δ x + ‖ B ‖ ∞ e 1 n Δ y + ‖ C 1 ‖ ∞ 4 e 1 n ( Δ y ) 2 + ‖ D ‖ ∞ e 1 n − | c | e 1 n + | b | ‖ A ‖ ∞ e ̅ 1 n + 1 Δ x + ‖ B ‖ ∞ e ̅ 1 n + 1 Δ y + ‖ C 1 ‖ ∞ 4 e ̅ 1 n + 1 ( Δ y ) 2 + ‖ D ‖ ∞ e ̅ 1 n + 1 − e ̅ 1 n + 1 .

Re-write inequality (35) as

(36) e 1 n + 1 ≤ 1 2 e Δ t + c ̅ 1 | a | + c ̅ 2 | a | + 4 d ̅ 1 | a | + | a | ‖ D ‖ ∞ − | c | | a | e 1 n + 1 2 + c ̅ 1 | b | + c ̅ 2 | b | + 4 d ̅ 1 | b | + | b | ‖ D ‖ ∞ − | b | e ̅ 1 n + 1 .

By using inequality (33) into inequality (36), we obtain

(37) e 1 n + 1 ≤ μ 1 e 1 n + μ 2 μ 3 e 1 n + Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) ,

where μ 1 = 1 2 e Δ t + c ̅ 1 | a | + c ̅ 2 | a | + 4 d ̅ 1 | a | + | a | ‖ D ‖ ∞ − | c | | a |

μ 2 = 1 2 + c ̅ 1 | b | + c ̅ 2 | b | + 4 d ̅ 1 | b | + | b | ‖ D ‖ ∞ − | b |

μ 3 = e Δ t + c ̅ 1 + c ̅ 2 + 4 d ̅ 1 + ‖ D ‖ ∞ ( e Δ t − 1 )

Let μ = μ 1 + μ 2 μ 3 then inequality (37) is expressed as

(38) e 1 n + 1 ≤ μ e 1 n + Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) ,

Put n = 0 in inequality (38) as

(39) e 1 1 ≤ μ e 1 n + Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) ,

Since e 1 0 = 0 due to the exact initial condition, so inequality (39) is written in the form

(40) e 1 1 ≤ Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) .

Substituting n = 1 in inequality (37) using (40) gives the following equation:

(41) e 1 2 ≤ μ e 1 n + Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) ≤ ( μ + 1 ) Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) .

If this continues for finite n , then the following inequality is obtained:

(42) e 1 n ≤ ( μ n − 1 + … + μ + 1 ) Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) = 1 − μ n 1 − μ Q ( O ( ( Δ t ) 2 , ( Δ x ) 2 , ( Δ y ) 2 ) ) .

For large n , i.e. n → ∞ the series … + μ n − 1 + … + μ + 1 becomes an infinite geometric series that will converge if ∣ μ ∣ < 1 .

4 Problem formulation

Examine the behaviour of a two-dimensional Williamson fluid that is laminar, incompressible, unstable, and flows over flat and oscillatory sheets. Let the x ⁎ -axis be placed horizontally, whereas the y ⁎ -axis is chosen to be perpendicular to the x ⁎ -axis. The streamwise coordinate of the flow is the x ⁎ -axis and cross-streamwise is the y ⁎ -axis. The flow in the fluid is generated by the sudden movement of the plate towards the positive direction of the x ⁎ -axis. The fluid is electrically conductive. The transverse magnetic field B ⇀ = ( 0 , B 0 , 0 ) and uniform electric field E ⇀ = ( 0 , 0 , − E 0 ) are applied in the flow region. It is to be noted that an electric field is stronger than a magnetic field, and a magnetic field follows Ohm’s law. J ˜ = σ ( E ⇀ + v ˜ × B ⇀ ) , where σ is the electric conductivity, J ˜ is a Joule current, and v ˜ is the velocity of fluid. The governing equation of this phenomenon can be written as

(43) ∂ u ⁎ ∂ x ⁎ + ∂ v ⁎ ∂ y ⁎ = 0 ,

(44) ∂ u ⁎ ∂ t ⁎ + u ⁎ ∂ u ⁎ ∂ x ⁎ + v ⁎ ∂ u ⁎ ∂ y ⁎ = ν ∂ 2 u ⁎ ∂ y ⁎ 2 + 2 ν Γ ∂ u ⁎ ∂ y ⁎ ∂ 2 u ⁎ ∂ y ⁎ 2 + σ ρ ( E 0 B 0 − B 0 2 u ⁎ ) ,

(45) ∂ T ∂ t ⁎ + u ⁎ ∂ T ∂ x ⁎ + v ⁎ ∂ T ∂ y ⁎ = α ∂ 2 T ∂ y ⁎2 + μ ρ C p ∂ u ⁎ ∂ y ⁎ 2 + μ ρ C p Γ ∂ u ⁎ ∂ y ⁎ 3 + σ ρ C p ( u ⁎ B 0 − E 0 ) 2 + 1 ρ C p q ′ ′ ′ ,

(46) ∂ C ∂ t ⁎ + u ⁎ ∂ C ∂ x ⁎ + v ⁎ ∂ C ∂ y ⁎ = D 1 ∂ 2 C ∂ y ⁎2 − k 1 ( C − C ∞ ) .

Subject to the boundary conditions

(47) u ⁎ = 0 , v ⁎ = 0 , T = 0 , C = 0 when t ⁎ = 0 , x ⁎ , y ⁎ > 0 u ⁎ = u w , v ⁎ = 0 , T = T w , C = C w when y ⁎ = 0 , x ⁎ , t ⁎ > 0 u ⁎ → 0 , T → T ∞ , C → C ∞ when y ⁎ → ∞ , x ⁎ , t ⁎ > 0 u ⁎ = 0 , v ⁎ = 0 , T = 0 , C = 0 when x ⁎ = 0 ,

where q ′ ′ ′ = k u w x ν ρ C p ( A ⁎ u ( T w − T ∞ ) + B ⁎ ( T − T ∞ ) ) represents the space- and temperature-dependent internal heat sources., A ⁎ and B ⁎ are spatial and temperature-dependent parameters, respectively. Γ is the time constant, ρ is the density of the fluid, σ represents the electrical conductivity, α is the thermal diffusivity, k 1 is the reaction rate, and D 1 is the mass diffusivity.

To reduce Eqs. (43)–(47) into dimensionless PDEs, consider the following transformations:

(48) u = u ⁎ u w , v = v ⁎ u w , t = u w t ⁎ L , y = y ⁎ L , x = x ⁎ L .

By employing transformations to Eqs. (43)–(47), the system of PDEs is given as follows:

(49) ∂ u ∂ x + ∂ v ∂ y = 0 ,

(50) ∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y = 1 Re ∂ 2 u ∂ y 2 + 2 We ∂ u ∂ y ∂ 2 u ∂ y 2 + M 2 ( E 1 − u ) ,

(51) ∂ θ ∂ t + u ∂ θ ∂ x + v ∂ θ ∂ y = 1 Pr 1 Re ∂ 2 θ ∂ y 2 + Ec Re ∂ u ∂ y 2 + EC W E ∂ u ∂ y 3 + EC M 2 ( u − E 1 ) 2 + L x 1 Pr ( A ⁎ u + B ⁎ θ ) ,

(52) ∂ ϕ ∂ t + u ∂ ϕ ∂ x + v ∂ ϕ ∂ y = 1 Sc 1 Re ∂ 2 ϕ ∂ y 2 − γ ϕ .

With the following dimensionless boundary conditions

(53) u = 0 , v = 0 , θ = 0 , ϕ = 0 when t = 0 , x , y > 0 u = 1 , v = 0 , θ = 1 , ϕ = 1 when y = 0 , x , t > 0 u → 0 , θ → 0 , ϕ → 0 when y → ∞ , x , t > 0 u = 0 , v = 0 , θ = 0 , ϕ = 0 when x = 0 ,

where We is the Weisenberg number, E 1 is the local electrical parameter, M is the magnetic parameter, Pr is the Prandtl number, Sc is the Schmidt number, γ is a dimensionless chemical reaction, Ec is the Exkert number, and Re is the Reynolds number, and these are expressed as follows:

We = ν Γ L 2 , E 1 = E 0 u w B 0 , M 2 = σ L B 0 2 ρ u w , Pr = ν α , Sc = ν D 1 , γ = L k 1 u w , Ec = u w 2 C p ( T w − T ∞ ) , Re = L u w ν

The following physical quantities indicated in the problem formulation are significant for better comprehension of the article:

Weisenberg number ( We ) : The relaxation time of the fluid as a function of the flow’s characteristic time scale is represented by this dimensionless parameter. It measures how much of an impact in non-Newtonian fluid flows is caused by viscoelastic vs inertial forces. The Weisenberg number increases as the impact of viscoelasticity on the flow behaviour becomes more significant.

Local electrical parameters ( E 1 ) : It helps determine the electric force level in any fluid as a part of the fluid flow process. Force on a fluid can be determined by knowing the different parameters related to the fluid. Another factor on which the movement of fluid depends extensively is the magnetic forces.

Magnetic parameter ( M 2 ) : In a fluid flow, this dimensionless quantity measures the inertial to electromagnetic forces ratio. Important for understanding the effect of magnetic fields on fluid dynamics, the parameter measures the relative strength of the field with respect to the fluid’s velocity.

Prandtl number ( Pr ) : This dimensionless variable represents the fluid’s momentum-to-thermal diffusivity ratio. It determines the thickness of the thermal boundary layer, flow heat transfer rates, and the importance of momentum vs thermal transport.

Schmidt number ( Sc ) : This metric quantifies the relationship between the diffusion of momentum and the diffusion of mass in the fluid. It measures the importance of moving momentum compared to moving mass and regulates the thickness of the layer where concentration changes occur and the speed at which mass is transferred in the flow.

Eckert number ( Ec ) : It is a measure of the correlation between the change in enthalpy and the kinetic energy of the fluid flow. The relative relevance of kinetic and thermal energy can be better understood by gaining insight into the flow’s equilibrium between kinetic energy dissipation and heat transfer.

The skin friction coefficients local Nusselt and Sherwood numbers are defined as follows:

(54) C f = τ w ρ u w 2 where τ w = μ ∂ u ⁎ ∂ y ⁎ + Γ 2 ∂ u ⁎ ∂ y ⁎ 2 y ⁎ = 0 Nu L = L q w k ( T w − T ∞ ) , q w = − k ∂ T ∂ y ⁎ y ⁎ = 0 Sh L = L q j D ( C w − C ∞ ) , q j = − D 1 ∂ C ∂ y ⁎ y ⁎ = 0 .

To make Eq. (54) dimensionless, transformations from Eq. (48) are substituted into Eq. (54), the dimensionless Skin friction coefficient, local Nusselt and Sherwood numbers are reduced as

(55) C f = 1 Re ∂ u ∂ y + We 2 ∂ u ∂ y 2 y = 0 Nu L = − ∂ θ ∂ y y = 0 Sh L = − ∂ ϕ ∂ y y = 0 .

5 Results and discussion

This contribution is based on the numerical scheme to discretize time-dependent problems. The scheme can be used to discretize those problems that contain only first-order partial derivatives in time. The scheme is a two-stage predictor–corrector stage. The first stage, or predictor stage, is first-order accurate in time, while the second stage, or corrector stage, provides second-order accuracy. The spatial discretization is performed by applying second-order central difference formulas. The scheme is conditionally stable, so the time step is restricted. This condition also depends on the choice of the parameter(s) in the considered problem(s). So, the scheme will remain stable if the appropriate numerical value(s) of the parameter(s) and suitable step sizes in space and time are chosen. The scheme is also consistent because it is constructed using the Taylor series and is second-order accurate. By Lax’s equivalence theorem, the scheme will converge for linear time-dependent problems. Since the continuity equation for incompressible flow does not contain the time derivative term, it is discretized by only first-order Euler’s scheme in both spatial coordinates. The proposed scheme is employed for Navier–Stokes’ energy and concentration equations.

Figure 1 illustrates the influence of the magnetic parameter on the velocity profile. The velocity profile decreases as the magnetic parameter increases. The drop in the velocity profile is caused by the rising Lorentz force that acts against the motion of the particles in the fluid. Figure 2 illustrates the impact of the Weisenberg number on the velocity profile. The velocity profile decreases as the Weisenberg number increases. An increase in the Weisenberg number results in higher fluid viscosity, leading to more excellent resistance to deformation. The presence of resistance causes a decrease in the flow speed, reducing the velocity profile. Figures 3 and 4 illustrate the impact of the coefficients of spatial and temperature-dependent factors on internal heat generation. The temperature profile increases as these values grow. Increasing these dimensionless factors intensifies the heat-creation process, resulting in an elevation in temperature. Figure 5 illustrates the impact of the Eckert number on the temperature distribution. The temperature profile becomes steeper as the Eckert number increases. With an increase in the Eckert number, both self-heating and the temperature profile rise due to the heightened dissipation resulting from the fluid’s internal friction. By applying the rules of electromagnetic and fluid dynamics, one may see why an electric field improves the velocity and temperature profile while a magnetic field has the opposite effect.

$Figure 1 Variation of magnetic parameter on velocity profile using Re = 3 , We = 0.01 , E 1 = 0.01 \text{Re}=3,\hspace{.25em}\text{We}=0.01,\hspace{.25em}{E}_{1}=0.01 .$

Figure 1

Variation of magnetic parameter on velocity profile using Re = 3 , We = 0.01 , E 1 = 0.01 .

$Figure 2 Variation of Weisenberg number on velocity profile using Re = 3 , M = 0.1 , E 1 = 0.01 \text{Re}=3,\hspace{.25em}M=0.1,\hspace{.25em}{E}_{1}=0.01 .$

Figure 2

Variation of Weisenberg number on velocity profile using Re = 3 , M = 0.1 , E 1 = 0.01 .

$Figure 3 Variation of coefficient of space-dependent term in a heat source for temperature profile using Re = 3 , We = 0.01 , M = 0.1 , \text{Re}=3,\hspace{.25em}\text{We}=0.01,\hspace{.25em}M=0.1,\hspace{.25em} E 1 = 0.01 , Pr = 0.9 , B ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 , γ = 0.1 {E}_{1}=0.01,\hspace{.25em}\text{Pr}=0.9,\hspace{.25em}{B}^{\ast }=0.1,\text{Ec}=0.1,\hspace{.25em}\text{Sc}=0.9,\hspace{.25em}\gamma =0.1 .$

Figure 3

Variation of coefficient of space-dependent term in a heat source for temperature profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , Pr = 0.9 , B ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 , γ = 0.1 .

$Figure 4 Variation of coefficient of the temperature-dependent term in a heat source for temperature profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , \text{Re}=3,\hspace{.25em}\text{We}=0.01,M=0.1,\hspace{.25em}{E}_{1}=0.01,\hspace{.25em} Pr = 0.9 , A ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 , γ = 0.1 \text{Pr}=0.9,\hspace{.25em}{A}^{\ast }=0.1,\text{Ec}=0.1,\hspace{.25em}\text{Sc}=0.9,\hspace{.25em}\gamma =0.1 .$

Figure 4

Variation of coefficient of the temperature-dependent term in a heat source for temperature profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , Pr = 0.9 , A ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 , γ = 0.1 .

$Figure 5 Variation of Eckert number on temperature profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , \text{Re}=3,\text{ We}=0.01,\text{ }M=0.1,\text{ }{E}_{1}=0.01, Pr = 0.9 , A ⁎ = 0.1 , B ⁎ = 0.1 , \text{Pr}=0.9,\text{ }{A}^{\text{⁎}}=0.1,{B}^{\text{⁎}}=0.1, Sc = 0.9 , \text{ Sc}=0.9, γ = 0.1 \gamma =0.1 .$

Figure 5

Variation of Eckert number on temperature profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , Pr = 0.9 , A ⁎ = 0.1 , B ⁎ = 0.1 , Sc = 0.9 , γ = 0.1 .

5.1 Adverse effect of magnetic field on velocity and temperature

5.1.1 Lorentz force

The Lorentz force is a force that acts on electrically conducting fluids as they pass through a magnetic field. The magnetic field and the current’s (or velocity’s) direction are perpendicular to this force. Consequently, it acts as a barrier to the fluid’s flow, causing the velocity to decrease.

5.1.2 Reduced energy transfer

The Lorentz force can also generate eddy currents within the fluid, leading to internal friction and the waste of kinetic energy. This dissipation reduces temperature, as energy is transformed into heat due to the resistance experienced by the fluid particles.

5.2 Enhancing effect of electric field

5.2.1 Electrothermal effects

In contrast to the magnetic field, the electric field can generate Joule heating in the fluid, resulting in a temperature rise. The heating effect occurs due to the fluid’s resistance to the passage of electric current.

5.2.2 Electrohydrodynamic effects

The electrohydrodynamic force describes the action of an electric field on charged particles in a fluid in specific situations. Under the right circumstances, this force can propel fluid motion while simultaneously amplifying the velocity field.

The effect of the Schmidt number on the concentration profile is shown in Figure 6. A higher Schmidt number is associated with a flatter concentration profile. The mass diffusivity drops, and the concentration profile flattens as the Schmidt number rises. Figure 7 shows how the concentration profile changes when the reaction rate parameter changes. When the reaction rate parameter is increased, the concentration profile gets flatter. Reactions that break or create chemical bonds between atoms transform concentration atoms into other substances, which causes the concentration profile to decay due to the process.

$Figure 6 Variation of Schmidt number on concentration profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , \text{Re}=3,\text{ We}=0.01,\text{ }M=0.1,\text{ }{E}_{1}=0.01, Pr = 0.9 , A ⁎ = 0.1 , \text{Pr}=0.9,\text{ }{A}^{\text{⁎}}=0.1, B ⁎ = 0.1 , {B}^{\text{⁎}}=0.1, Ec = 0.1 , γ = 0.1 \text{ Ec}=0.1,\text{ }\gamma =0.1 .$

Figure 6

Variation of Schmidt number on concentration profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , Pr = 0.9 , A ⁎ = 0.1 , B ⁎ = 0.1 , Ec = 0.1 , γ = 0.1 .

$Figure 7 Variation of reaction rate parameter on concentration profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , \text{Re}=3,\text{ We}=0.01,\text{ }M=0.1,\text{ }{E}_{1}=0.01, Pr = 0.9 , A ⁎ = 0.1 , \text{Pr}=0.9,\text{ }{A}^{\text{⁎}}=0.1, B ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 {B}^{\text{⁎}}=0.1,\text{ Ec}=0.1,\text{ Sc}=0.9 .$

Figure 7

Variation of reaction rate parameter on concentration profile using Re = 3 , We = 0.01 , M = 0.1 , E 1 = 0.01 , Pr = 0.9 , A ⁎ = 0.1 , B ⁎ = 0.1 , Ec = 0.1 , Sc = 0.9 .

Figure 8 illustrates the impact of the magnetic parameter and Weisenberg number on the skin friction coefficient. The skin friction coefficient decreases as the magnetic parameter and Weisenberg number increase. This decline in skin friction coefficient results from a reduced velocity profile due to magnetic parameters and Weisenberg number growth. So, the skin friction rises when the velocity of the flow grows. Figure 9 displays the effect of Prandtl and Eckert numbers on local Nusselt numbers. Skin friction coefficient rises and decays by enhancing the Prandtl and Eckert numbers, respectively. Since heat transfer is better for growing thermal conductivity, conduction heat transfer is affected by the declining thermal conductivity that enhances the local Nusselt number. Figure 10 illustrates the impact of the Schmidt number and reaction rate parameter on the local Sherwood number. The local Sherwood number rises as the Schmidt number and reaction rate parameter increment. Since mass diffusivity decays by rising Schmidt number, this leads to a decline in diffusion rate, and so local Sherwood number escalates. Figures 11–14 show the contour plots of horizontal and vertical components of the velocity profile over space and time coordinates for the oscillatory boundary conditions. Two different types of oscillatory conditions are considered with different periods. Its effect along the time coordinates on the horizontal component of velocity can be seen in Figures 11–14.

$Figure 8 Variation of magnetic parameter and Weisenberg number on skin friction coefficient using Re = 1 , E 1 = 0.01 \text{Re}=1,\hspace{.25em}{E}_{1}=0.01 .$

Figure 8

Variation of magnetic parameter and Weisenberg number on skin friction coefficient using Re = 1 , E 1 = 0.01 .

$Figure 9 Variation of Prandtl number and Eckert number on local Nusselt number using Re = 1 , E 1 = 0.01 , We = 0.1 , M = 0.1 , A ⁎ = 0.1 , B ⁎ = 0.1 \text{Re}=1,\hspace{.25em}{E}_{1}=0.01,\hspace{.25em}\text{We}=0.1,\hspace{.25em}M=0.1,\hspace{.25em}{A}^{\ast }=0.1,\hspace{.25em}{B}^{\ast }=0.1 .$

Figure 9

Variation of Prandtl number and Eckert number on local Nusselt number using Re = 1 , E 1 = 0.01 , We = 0.1 , M = 0.1 , A ⁎ = 0.1 , B ⁎ = 0.1 .

$Figure 10 Variation of Schmidt number and reaction rate parameter on local Sherwood number using Re = 1 , E 1 = 0.01 , We = 0.1 , M = 0.1 \text{Re}=1,\hspace{.25em}{E}_{1}=0.01,\hspace{.25em}\text{We}=0.1,\hspace{.25em}M=0.1 .$

Figure 10

Variation of Schmidt number and reaction rate parameter on local Sherwood number using Re = 1 , E 1 = 0.01 , We = 0.1 , M = 0.1 .

$Figure 11 Contour plot for the horizontal component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , \text{Re}=1,\text{ }{E}_{1}=0.01,\text{ We}=0.01,\text{ }M=0.1,\text{ }{L}_{y}=25,\text{ }{L}_{x}=27, U w = 0.5 cos ( 1.5 t ) {U}_{\text{w}}=0.5\text{ cos}(1.5t) .$

Figure 11

Contour plot for the horizontal component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , U w = 0.5 cos ( 1.5 t ) .

$Figure 12 Contour plot for the vertical component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , \text{Re}=1,\text{ }{E}_{1}=0.01,\text{ We}=0.01,\text{ }M=0.1,\text{ }{L}_{y}=25,\text{ }{L}_{x}=27, U w = 0.5 cos ( 1.5 t ) {U}_{\text{w}}=0.5\text{ cos}(1.5t) .$

Figure 12

Contour plot for the vertical component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , U w = 0.5 cos ( 1.5 t ) .

$Figure 13 Contour plot for the horizontal component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , \text{Re}=1,\text{ }{E}_{1}=0.01,\text{ We}=0.01,\text{ }M=0.1,\text{ }{L}_{y}=25,\text{ }{L}_{x}=27, U w = 0.5 ( cos ( 1.5 t ) + sin ( 1.5 t ) ) {U}_{\text{w}}=0.5(\text{cos}(1.5t)+\text{sin}(1.5t)) .$

Figure 13

Contour plot for the horizontal component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , U w = 0.5 ( cos ( 1.5 t ) + sin ( 1.5 t ) ) .

$Figure 14 Contour plot for the vertical component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , \text{Re}=1,\text{ }{E}_{1}=0.01,\text{ We}=0.01,\text{ }M=0.1,\text{ }{L}_{y}=25,\text{ }{L}_{x}=27, U w = 0.5 ( cos ( 1.5 t ) + sin ( 1.5 t ) ) {U}_{\text{w}}=0.5(\text{cos}(1.5t)+\text{sin}(1.5t)) .$

Figure 14

Contour plot for the vertical component of velocity profile using Re = 1 , E 1 = 0.01 , We = 0.01 , M = 0.1 , L y = 25 , L x = 27 , U w = 0.5 ( cos ( 1.5 t ) + sin ( 1.5 t ) ) .

Table 1 compares the suggested scheme and the second-order Runge–Kutta method. The comparison is based on the norm of error and the time consumption, for Example 1, as discussed in Li et al. [37]. The proposed scheme gives a smaller error compared with the second-order Runge–Kutta method, but the chosen Runge–Kutta method is computationally less expensive than the proposed scheme.

Table 1

Comparison of proposed scheme with existing scheme using N x = 50 = N y ( no . grid points ) , t f ( final time ) = 0.1

Δt	Proposed L 2 error	Time (s)	Runge–Kutta L 2 error	Time (s)
0.1 / 1 , 500	2.7479 × 10⁻⁴	28.9	5.3073 × 10⁻⁴	27.6
0.1 / 1 , 750	2.7472 × 10⁻⁴	34.2	5.0873 × 10⁻⁴	32.6
0.1 / 2 , 000	2.7470 × 10⁻⁴	38.9	4.9222 × 10⁻⁴	36.3
0.1 / 2 , 250	2.7467 × 10⁻⁴	42.3	4.7939 × 10⁻⁴	40.9

6 Conclusion

This study introduced a new method for efficiently solving time-dependent PDEs using an exponential integrator. The stability and second-order precision of the scheme in time have been proven by Fourier series analysis, proving its durability as a numerical tool for solving complicated fluid dynamics problems. We have learned a lot about the system’s behaviour by applying this technique to model and simulate the flow of Williamson fluid under the effect of a space and temperature-dependent heat source. A complete framework for studying the complex interaction of parameters has been developed by discretizing the continuity equation using a first-order scheme and applying the suggested exponential integrator to the Navier–Stokes, energy, and concentration equations with central difference in space. The complex consequences of changing the coefficients of the space and temperature terms of the heat source under consideration have been better understood by the display of concentration, velocity, and temperature profiles under various situations. A proposed computational scheme and its application on non-Newtonian fluid have been provided. The scheme consisted of two stages. The first stage is the first-order exponential integrator, and the second stage corrects the solution. Both stages combined proved the second-order exponential solution accurate. The effects of these variables on the Sherwood number, local Nusselt number, and skin friction coefficient are shown graphically. The results of the investigation led to the following findings:

A decrease in the velocity profile was seen as the Weisenberg number and magnetic parameter increased.
The increasing magnetic parameter and Weisenberg number caused a decrease in the skin friction coefficient.
The temperature profile has increased by increasing the coefficients of the spatial and temperature-dependent factors of the heat source.

The proposed exponential integrator also captures the system’s behavior, making it a potentially valuable numerical tool for investigating and understanding such complicated processes. The scheme’s explicit character and second-order precision make it a good choice for modelling time-dependent PDEs in many scientific and engineering applications.

This study concludes that new numerical approaches and sophisticated physical system investigations can help us understand complex fluid dynamics [38–40]. This discovery advances non-Newtonian fluid flow theory and makes the exponential integrator more feasible for dynamic and heterogeneous heat source situations. This work lays the framework for future numerical method developments and their application to more complex physical processes, with implications beyond the system explored here.

Acknowledgments

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

Funding information: The authors state no funding involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Williamson RV. The flow of pseudoplastic materials. Ind Eng Chem. 1929;21:1108–11.10.1021/ie50239a035Suche in Google Scholar

[2] Nadeem S, Hussain ST. Heat transfer analysis of williamson fluid over exponentially stretching surface. Appl Math Mech. 2014;35:489–502.10.1007/s10483-014-1807-6Suche in Google Scholar

[3] Amjad M, Ahmed K, Akbar T, Muhammad T, Ahmed I, Alshomrani AS. Numerical investigation of double diffusion heat flux model in williamson nanofluid over an exponentially stretching surface with variable thermal conductivity. Case Stud Therm Eng. 2022;36:102231.10.1016/j.csite.2022.102231Suche in Google Scholar

[4] Kothandapani M, Prakash J. Effects of thermal radiation parameter and magnetic field on the peristaltic motion of williamson nanofluids in a tapered asymmetric channel. Int J Heat Mass Transf. 2015;81:234–45.10.1016/j.ijheatmasstransfer.2014.09.062Suche in Google Scholar

[5] Ahmed K, Akbar T. Numerical investigation of magnetohydrodynamics williamson nanofluid flow over an exponentially stretching surface. Adv Mech Eng. 2021;13:16878140211019876.10.1177/16878140211019875Suche in Google Scholar

[6] Hayat T, Bashir G, Waqas M, Alsaedi A. MHD 2D flow of williamson nanofluid over a nonlinear variable thicked surface with melting heat transfer. J Mol Liq. 2016;223:836–44.10.1016/j.molliq.2016.08.104Suche in Google Scholar

[7] Abbas A, Khan A, Abdeljawad T, Aslam M. Numerical simulation of variable density and magnetohydrodynamics effects on heat generating and dissipating Williamson Sakiadis flow in a porous space: Impact of solar radiation and Joule heating. Heliyon. 2023;9(11):e21726.10.1016/j.heliyon.2023.e21726Suche in Google Scholar PubMed PubMed Central

[8] Shafiq A, Çolak AB, Sindhu TN, Al-Mdallal QM, Abdeljawad T. Estimation of unsteady hydromagnetic Williamson fluid flow in a radiative surface through numerical and artificial neural network modeling. Sci Rep. 2021;11(1):14509.10.1038/s41598-021-93790-9Suche in Google Scholar PubMed PubMed Central

[9] Ahmed K, Akbar T, Muhammad T. Physical aspects of homogeneous-heterogeneous reactions on MHD williamson fluid flow across a nonlinear stretching curved surface together with convective boundary conditions. Math Probl Eng. 2021;2021:7016961.10.1155/2021/7016961Suche in Google Scholar

[10] Ahmed K, Khan WA, Akbar T, Rasool G, Alharbi SO, Khan I. Numerical investigation of mixed convective williamson fluid flow over an exponentially stretching permeable curved surface. Fluids. 2021;6:260.10.3390/fluids6070260Suche in Google Scholar

[11] Makinde OD, Mabood F, Ibrahim MS. Chemically reacting on MHD boundary-layer flow of nanofluids over a nonlinear stretching sheet with heat source/sink and thermal radiation. Therm Sci. 2018;22:495–506.10.2298/TSCI151003284MSuche in Google Scholar

[12] Turkyilmazoglu M. An analytical treatment for the exact solutions of MHD flow and heat over two--three dimensional deforming bodies. Int J Heat Mass Transf. 2015;90:781–9.10.1016/j.ijheatmasstransfer.2015.07.025Suche in Google Scholar

[13] Mabood F, Ibrahim SM, Lorenzini G. Chemical reaction effects on MHD rotating fluid over a vertical plate embedded in porous medium with heat source. J Eng Thermophys. 2017;26:399–415.10.1134/S1810232817030109Suche in Google Scholar

[14] Hashemi-Tilehnoee M, Dogonchi AS, Seyyedi SM, Chamkha AJ, Ganji DD. Magnetohydrodynamic natural convection and entropy generation analyses inside a nanofluid-filled incinerator-shaped porous cavity with wavy heater block. J Therm Anal Calorim. 2020;141:2033–45.10.1007/s10973-019-09220-6Suche in Google Scholar

[15] Rashidi S, Esfahani JA, Maskaniyan M. Applications of magnetohydrodynamics in biological systems – review on the numerical studies. J Magn Magn Mater. 2017;439:358–72.10.1016/j.jmmm.2017.05.014Suche in Google Scholar

[16] Madhu M, Shashikumar NS, Thriveni K, Gireesha BJ, Mahanthesh B. Irreversibility analysis of the MHD williamson fluid flow through a microchannel with thermal radiation. Waves Random Complex Media. 2022;1–23. 10.1080/17455030.2022.2111473.Suche in Google Scholar

[17] Mishra P, Kumar D, Kumar J, Abdel-Aty AH, Park C, Yahia IS. Analysis of MHD Williamson Micropolar Fluid Flow in Non-Darcian Porous Media with Variable Thermal Conductivity. Case Stud. Therm Eng. 2022;36:102195.10.1016/j.csite.2022.102195Suche in Google Scholar

[18] Almaneea A. Numerical study on heat and mass transport enhancement in MHD williamson fluid via hybrid nanoparticles. Alex Eng J. 2022;61:8343–54.10.1016/j.aej.2022.01.041Suche in Google Scholar

[19] Reddy MV, Lakshminarayana P. MHD radiative flow of williamson nanofluid with cattaneo-christov model over a stretching sheet through a porous medium in the presence of chemical reaction and suction/injection. J Porous Media. 2022;25:1–15.10.1615/JPorMedia.2022041423Suche in Google Scholar

[20] Asjad MI, Zahid M, Inc M, Baleanu D, Almohsen B. Impact of activation energy and MHD on williamson fluid flow in the presence of bioconvection. Alex Eng J. 2022;61:8715–27.10.1016/j.aej.2022.02.013Suche in Google Scholar

[21] Raju CSK, Sandeep N, Babu MJ, Sugunamma V. Dual solutions for three-dimensional MHD flow of a nanofluid over a nonlinearly permeable stretching sheet. Alex Eng J. 2016;55(1):151–62. 10.1016/j.aej.2015.12.017 Suche in Google Scholar

[22] Mather DJ. The motion of viscous liquid past a paraboloid. Quart J Mech Appl Math. 1961;14(4):423–30. 10.1093/qjmam/14.4.423 Suche in Google Scholar

[23] Lee LL. Boundary layer over a thin needle. Phys Fluids. 1967;10:820. 10.1063/1.1762194 Suche in Google Scholar

[24] Cebeci T, Na TY, Mosinskis G. Laminar boundary layers on slender paraboloids. AIAA J. 1969;7(1):1372–4.10.2514/3.5353Suche in Google Scholar

[25] Miller DR. The boundary layer on a paraboloid of revolution. Math Proc Camb Philos Soc. 1969;65(1):285–99.10.1017/S0305004100044248Suche in Google Scholar

[26] Miller DR. The downstream solution for steady viscous flow past a paraboloid. Math Proc Camb Philos Soc. 1971;70(1):123–33. 10.1017/s0305004100049732 Suche in Google Scholar

[27] Veldman AEP. The numerical solution of the Navier-Stokes equations for laminar incompressible flow past a paraboloid of revolution. Comput Fluids. 1973;1(3):251–71. 10.1016/0045-7930(73)90010-8 Suche in Google Scholar

[28] Nadeem S, Ashiq S, Ali M. Williamson fluid model for the peristaltic flow of chyme in small intestine. Math Probl Eng. 2012;12(12):1–4.10.1155/2012/479087Suche in Google Scholar

[29] Shehzad SA, Abdullah Z, Alsaedi A, Abbasi FM, Hayat T. Thermally radiative three-dimensional flow of Jeffrey nanofluid with internal heat generation and magnetic field. J Magn Magn Mater. 2015;397:108–14. 10.1016/j.jmmm.2015.07.057 Suche in Google Scholar

[30] Motsa SS, Animasaun IL. Paired quasi-linearization analysis of heat transfer in unsteady mixed convection nanofluid containing both nanoparticles and gyrotactic microorganisms due to impulsive motion. J Heat Transf. 2016;138(11):114503. 10.1115/1.4034039 Suche in Google Scholar

[31] Abbasi FM, Shehzad SA, Hayat T, Ahmad B. Doubly stratified mixed convection flow of Maxwell nanofluid with heat generation/absorption. J Magn Magn Mater. 2015;404:159–65. 10.1016/j.jmmm.2015.11.090 Suche in Google Scholar

[32] Sandeep N, Koriko OK, Animasaun IL. Modified kinematic viscosity model for 3D-Casson fluid flow within boundary layer formed on a surface at absolute zero. J Mol Liq. 2016;221:1197–206. 10.1016/j.molliq.2016.06.049 Suche in Google Scholar

[33] Bhatti MM, Öztop HF, Ellahi R, Sarris IE, Doranehgard MH. Insight into the investigation of diamond (C) and Silica (SiO2) nanoparticles suspended in water-based hybrid nanofluid with application in solar collector. J Mol Liq. 2022;357:119134.10.1016/j.molliq.2022.119134Suche in Google Scholar

[34] Bhatti MM, Shahid A, Sarris IE, Anwar Bég O. Spectral relaxation computation of Maxwell fluid flow from a stretching surface with quadratic convection and non-Fourier heat flux using Lie symmetry transformations. Int J Mod Phys B. 2023;37(09):2350082.10.1142/S0217979223500820Suche in Google Scholar

[35] Nadeem S, Ishtiaq B, Alzabut J, Ghazwani HA. Entropy generation for exact irreversibility analysis in the MHD channel flow of Williamson fluid with combined convective-radiative boundary conditions. Heliyon. 202410.1016/j.heliyon.2024.e26432Suche in Google Scholar PubMed PubMed Central

[36] Nawaz Y, Arif MS, Abodayeh K, Soori AH, Javed U. A modification of explicit time integrator scheme for unsteady power-law nanofluid flow over the moving sheets. Front Energy Res. 2024;12:1335642.10.3389/fenrg.2024.1335642Suche in Google Scholar

[37] Li FL, Wu ZK, Ye CR. A finite difference solution to a two-dimensional parabolic inverse problem. Appl Math Model. 2012;36(5):2303–13.10.1016/j.apm.2011.08.025Suche in Google Scholar

[38] Arif MS, Abodayeh K, Nawaz Y. A finite difference explicit-implicit scheme for fractal heat and mass transportation of Williamson nanofluid flow in quantum calculus. Numer Heat Transf Part A. 2024;1–23. 10.1080/10407782.2024.2308753.Suche in Google Scholar

[39] Arif MS, Shatanawi W, Nawaz Y. Finite element study of electrical MHD Williamson nanofluid flow under the effects of frictional heating in the view of viscous dissipation. Energies. 2023;16(6):2778.10.3390/en16062778Suche in Google Scholar

[40] Arif MS, Abodayeh K, Nawaz Y. Innovative stochastic finite difference approach for modelling unsteady non-Newtonian mixed convective fluid flow with variable thermal conductivity and mass diffusivity. Front Phys. 2024;12:1373111.10.3389/fphy.2024.1373111Suche in Google Scholar

Received: 2024-01-20

Revised: 2024-04-20

Accepted: 2024-05-02

Published Online: 2024-05-30

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/phys-2024-0034

Schlagwörter für diesen Artikel

exponential integrator; stability; convergence; Williamson fluid; heat source

Creative Commons

BY 4.0