A high-accuracy exponential time integration scheme for the Darcy–Forchheimer Williamson fluid flow with temperature-dependent conductivity

Muhammad Shoaib Arif

doi:10.1515/phys-2025-0190

Article Open Access

A high-accuracy exponential time integration scheme for the Darcy–Forchheimer Williamson fluid flow with temperature-dependent conductivity

Muhammad Shoaib Arif

Published/Copyright: July 25, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 23 Issue 1

Abstract

Accurately modelling non-Newtonian fluid flow through porous media is vital in several industrial and engineering processes, particularly where inertial effects and temperature-dependent properties are critical. This study aims to develop a high-order numerical scheme for simulating Darcy–Forchheimer flow of Williamson fluids with temperature-dependent thermal conductivity. The proposed methodology combines a second-order exponential time integrator with a sixth-order compact finite difference scheme for spatial discretization. A modified predictor–corrector structure is formulated, and its conditional stability is verified through Fourier (von Neumann) analysis. The method avoids linearization and iterative solvers, enhancing computational efficiency. Quantitative validation demonstrates that the scheme achieves second-order accuracy in time and sixth-order convergence in space, with a reduction in error norm of up to 18% compared to the classical second-order Runge–Kutta method. Parametric studies confirm the strong influence of Weissenberg number, Forchheimer number, and non-linear conductivity on velocity, temperature, and concentration fields, validating the robustness and applicability of the scheme in complex transport phenomena. All computations were carried out using MATLAB R2023a, where the proposed scheme was implemented and validated against benchmark solutions to confirm its numerical accuracy and computational efficiency.

Keywords: computational scheme; stability analysis; convergence behaviour; Darcy–Forchheimer flow; Williamson fluid; heat and mass transfer; porous media; mixed convection; numerical comparison

Nomenclature

u: velocity component in the x -direction ( m / s )
y: spatial coordinate perpendicular to the plate ( m )
T: temperature of the fluid ( K )
C: concentration ( kg / m³ )
t: time ( s )
ν: kinematic viscosity ( m² / s )
ρ: density ( kg / m³ )
μ: dynamic viscosity ( Pa·s ) a
B ∘: strength of the magnetic field ( T (tesla))
k p: permeability of porous medium ( m 2 )
σ: electrical conductivity ( S / m )
c p: specific heat capacity ( J / ( kg·K ) )
k ∞: thermal conductivity ( W / ( m·K ) )
Γ: fluid relaxation time (Williamson) ( s )
g: gravitational acceleration ( m / s² )
β T: thermal expansion coefficient ( 1 / K )
β c: solutal expansion coefficient ( 1 / ( kg / m ³ ) )
D: mass diffusivity ( m² / s )
k 1: dimensional reaction rate parameter ( 1 / s )
b ̅: forchheimer inertial coefficient ( 1 / m )
A *: space-dependent heat generation parameter
B *: temperature-dependent heat generation parameter

Dimensionless quantities

θ: temperature
ϕ: concentration
W e: Weissenberg number
D a: Darcy number
F s: Forchheimer number
M: magnetic field parameter
P r: Prandtl number
S c: Schmidt number
E c: Eckert number
N: buoyancy ratio parameter
K c: reaction rate
ϵ 1: thermal conductivity variation parameter
ϵ 2: amplitude parameter for boundary oscillations

Subscripts

∞: ambient/reference condition
w: wall condition
L: local (e.g. local Nusselt or Sherwood number)

Superscripts

*: dimensional variable

1 Introduction

The fluid dynamics investigation within porous media has attracted considerable interest due to its extensive applicability in engineering and industrial processes. Darcy–Forchheimer flow, which enhances the traditional Darcy law by including inertial factors, is essential for explaining complicated fluid transport phenomena. The complexity of mathematical modelling increases when integrated with non-Newtonian fluid behaviour, exemplified by Williamson fluids. These fluids demonstrate shear-thinning characteristics, making them suitable for modelling practical situations, including increased oil recovery, chemical processes, and biological applications.

Non-Newtonian fluids have many uses but encounter complex boundary conditions and interactions with porous materials. One of these applications is increased oil recovery, and another is geothermal systems, where porous media flow is an essential component. Ideally suited to these cases is the Darcy–Forchheimer model, which extends Darcy’s law and uses inertial effects to depict flows at high velocities in porous media. This model is constructive for characterizing flows in porous media since many applications involving high-velocity flows have a non-linear pressure–velocity link.

The impact of temperature-dependent thermal conductivity must be carefully considered when researching flows of non-Newtonian fluids through porous media. Systems with large temperature gradients, such as cooling systems, geothermal reservoirs, and polymer synthesis, necessitate temperature-dependent changes in thermal conductivity for accurate thermal transport modelling. The system’s actual temperature response must be considered. Hence, it is not always correct to assume constant thermal conductivity. By incorporating temperature-dependent thermal conductivity, the heat transport of the fluid can be more accurately and realistically represented.

The motivation of the study is the desire to precisely model complicated transport phenomena in porous and thermally active settings drives research on the Darcy–Forchheimer flow of Williamson fluids with temperature-dependent thermal conductivity. In many real-world applications, such as enhanced oil recovery, polymer extrusion, biomedical flows, and geothermal systems, fluids show non-Newtonian behaviour, interact with porous media, and react to temperature and solutal gradients under electromagnetic fields. Classical models and low-order numerical techniques often lack the necessary precision or computational efficiency to represent these coupled effects. Thus, this work intends to create a high-accuracy, two-stage numerical method that manages the governing equations’ natural non-linearities without using iterative solvers. The development and study of this model handle the flow system’s physical complexity and the mathematical requirement for stable, accurate, and efficient simulation tools.

1.1 Industrial significance and research gap

In many engineering and industrial applications, the mathematical modelling of non-Newtonian fluid flow through porous media is vital. Williamson fluids, especially those showing shear-thinning characteristics, are commonly used in increased oil recovery, polymer extrusion, biomedical transport (e.g. drug delivery and blood flow in tissues), and thermal management systems including nanofluid-based porous heat exchangers. Especially when temperature- and concentration-dependent impacts are present, these processes can include complex interactions involving viscous, inertial, magnetic, and buoyant forces. Designing effective and dependable systems in geothermal energy extraction, chemical reactors, and industrial filtration units depends on precisely capturing these influences.

Although many studies on Newtonian and certain non-Newtonian models have been published, little attention has been paid to the combined impacts of Williamson rheology, Darcy–Forchheimer porous drag, variable thermal conductivity, and transient boundary conditions. Furthermore, current numerical methods sometimes lack the necessary accuracy and computing efficiency to resolve such strongly coupled, non-linear systems. The current work offers a high-accuracy numerical framework combining a two-stage time-stepping approach with a sixth-order compact spatial discretization to close this gap. The model provides a strong basis for simulating real-world transport processes in porous domains, intended to handle mixed convective flow, non-linear thermal effects, and reactive mass transport.

The Williamson fluid is a prominent example of a pseudoplastic fluid model. The Williamson fluid model has been the subject of many published studies over the past 20 years. The rheological approaches employed with Williamson fluids to consider how the Weissenberg number affects the pumping properties and flow have been attempted by numerous scholars in these works. The Williamson fluid flow towards stretchable surfaces has captured the curiosity of multiple researchers in numerous vital sectors, including polymer extrusion, plastic films, metal spinning, and metallurgical operations, among many more. The outcomes of a model that Williamson [1] created in 1929 to assign the flow of pseudoplastic liquids were shown. Next, Nadeem and Akram [2] found the analytical solution for the peristaltic flow of a Williamson fluid in an asymmetric channel after showing that the pressure gradient reduces with increasing Weissenberg number. An analytical solution for the boundary layer flow of the Williamson fluid model approaching a stretched surface was provided in a separate publication by Nadeem et al. [3]. Investigating the peristaltic flow of a Williamson fluid in an asymmetric channel, Akbar et al. [4] employed numerical approaches. Their results show that the nanoparticle concentration decreases as the Brownian motion parameter increases. Considering the impacts of viscous dissipation and Joule heating, Eldabe et al. [5] examined the peristaltic flow of a Williamson fluid that conducts electricity through a porous medium. Their research shows that the concentration profile flattens out with increasing values of the magnetic field parameter. In the presence of a chemical reaction, Krishnamurthy et al. [6] numerically solved the Williamson nanofluid flow and melting heat transfer past a horizontal plate. They demonstrated that the solutal boundary barrier becomes thinner as the chemical reaction parameter increases.

In addition, Bhatti and Rashidi [7] discussed how thermo-diffusion and heat radiation affect Williamson nanofluid. According to their findings, increasing the Lewis number lowers the concentration profile. In their discussion of the impacts of a non-linear variable thickness surface on the maximum hydrodynamic drag (MHD) flow of Williamson fluid, Hayat et al. [8] focused on melting. They showed that the heat transfer rate increases as the thermophoresis parameter values increase. The recent investigation by Kumaran and Sandeep [9] investigates the parabolic flow characteristics of Williamson nanofluid, particularly highlighting the implications of cross-diffusion effects. Hayat et al. [10] used modified Darcy's law to conduct numerical investigations of Williamson fluids with an endoscope. Axial velocity was discovered to be reduced for larger Hartmann numbers. Not long ago, Hamid et al. investigated the unstable Williamson fluid flow caused by a wedge-shaped geometry [11]. According to their proposal, a larger wedge angle parameter results in a greater surface shear stress. The effects of wedge-geometrically driven Williamson fluid flow were studied by Hashim et al. [12], who showed that a greater unsteady parameter results in a faster heat transfer rate. Using a magnetic field and varying thermal conductivity, Hashim et al. [13] studied the Williamson nanofluid mixed convection flow across a radially extended surface. They successfully showed that as the thermal conductivity parameter increases, the fluid temperature and the thickness of the thermal boundary layer both increase. The energy efficiency and fluid flow mechanism of a non-Newtonian liquid with nanoparticles were shown by Bahiraei et al. [14]. Makinde et al. [15] used analytical methods to evaluate the effects of chemical processes and heat generation/absorption on the behaviour of an electrically conducting nanofluid affected by a non-linear stretched surface. Nanofluid heat transmission in a rotating system was quantitatively studied by Mabood et al. using thermophoresis and Brownian motion features [16].

In the presence of the Darcy–Forchheimer flow, the stratification and velocity slip conditions across the surface of the stretchable cylinder have been investigated [17]. The two-phase magneto-hydrodynamics flow in an entropy-optimized nanofluid via a stretched sheet in the presence of microorganisms and an extended Darcy–Forchheimer porous medium is studied analytically and numerically in this study [18]. Research on fluid dynamics, specific fluid structures, and their interactions with homogeneous porous media will be aided by the current results on thermally magnetized Williamson fluid flow fields obtained through symmetry analysis [19].

In chemical reactions, there are two types of interactions: homogeneous and heterogeneous. The catalyst concentration near the surface can be increased, and the bulk liquid concentration can be decreased by adjusting the heterogeneous reaction parameters. Several scholars have now proposed a definition of stretchable flows that considers both homogeneous and heterogeneous effects. In his initial study, Merkin [20] focused on the homogeneous–heterogeneous reactions in viscous fluid laminar flow on a flat plate. He found heterogeneous reactions on catalyst surfaces and homogeneous reactions in cubic autocatalysis. It is dominating close to the flat plate reaction’s leading edge.

Chaudhary and Merkin [21] introduced homogeneous–heterogeneous reactions exhibiting equal diffusivities. Chemically reactive species were studied by Ziabakhsh et al. [22] as they moved and diffused across a non-linear stretching surface in a porous medium. The results of homogeneous–heterogeneous reactions in a viscoelastic fluid moving across a stretched surface were reviewed by Khan and Pop [23]. The significant viscoelastic parameter diminished the surface concentration, according to their investigation. Kameswaran et al. [24] examined how a nanofluid’s porous stretching surface was affected by homogeneous–heterogeneous processes. They discovered that a weaker heterogeneous response weakens the concentration profile. The properties of carbon nanotubes in stagnation point flows caused by non-linear stretched surfaces of varying thicknesses were recently detailed by Hayat et al. [25]. According to their findings, the drag surface force reduces as the ratio parameter increases. Chen et al. [26] computationally resolved the homogeneous–heterogeneous interaction on the homogeneous ignition characteristic in hydrogen-fuelled catalytic microreactors operating in the submillimeter to millimetre range. Recent numerical discussion by Hashim and Khan [27] focused on Carreau fluids with homogeneous–heterogeneous processes. Joule heating and homogeneous–heterogeneous reactions pertaining to flow across a non-linear stretched surface were addressed by Khan et al. [28]. Furthermore, the characteristics of homogeneous–heterogeneous processes in the three-dimensional Sisko fluid flow over a bidirectional extending surface were examined by Khan et al. [29]. To discover a precise series solution for fluid flow problems, the study effort [30,31,32,33,34] presents an effective model based on the Levenberg–Marquardt algorithm and artificial neural networks (LMA-BANN).

Several models have explained the behaviour of non-Newtonian fluids. The Williamson model is one of the models included in the subclass. This model describes the dynamics of the shear-thinning characteristics, a concept developed by Williamson [35]. The time-dependent two-dimensional WF flow on a wedge was elucidated by Hashim et al. [36]. They illustrated that with an increase in the Weissenberg number, there is a corresponding decrease in the thermal state of the fluid. Khan et al. [37] examined the effects of alterations in the viscosity of MHD WNF on a non-linear SS. The findings indicate that the fluid velocity diminishes with an increase in the Weissenberg number. Salahuddin [38] analysed heat and mass transfer concerning WF on an oscillatory surface. He observed that the Weissenberg number results in a diminishing characteristic in the movement of the fluid particles. Bilal et al. [39] conducted a numerical investigation of the magnetohydrodynamic wall-free flow around a cylinder. It was observed that the FT increases with an increase in the Weissenberg number. The phenomenon of MHD WNF on a wedge was documented by Subbarayudu et al. [40]. They observed that the fluid diminishes its warmth due to the influence of the Weissenberg number.

In recent decades, numerous researchers have shown a keen interest in examining the non-linear thermal radiation effect, as it plays a crucial role in various industrial processes. Examples include polymer processing, missile technology, propulsion systems, satellites, solar energy innovations, and aircraft design. Given these considerations, Ganesh Kumar et al. [41] employed a stretched sheet methodology to investigate the three-dimensional flow of Jeffrey nanofluids. As the values of the radiation parameter are increased, there is a corresponding reduction in the HTG. Usman et al. [42] explored the effects of non-linear radiative Eyring Powell NF flow on a rotating plate influenced by gyrotactic microorganisms. The investigation revealed that the surface heating parameter contributes to the enhancement of the HTG. The complicated three-dimensional Jeffrey nanofluid non-linear radiative flow was explained by Shehzad et al. [43]. According to the observation, an increase in the radiation parameter is correlated with an increase in nanoparticle concentration. Hayat et al. [44] thoroughly analysed the consequences of three-dimensional nanofluids’ non-linear thermal radiative flow. The data show that the LNN improves as the radiation parameter increases. Usman et al. [45] studied micropolar nanofluids with activation energy in terms of their non-linear radiative flow. The effects of non-linear heat radiation on viscoelastic fluids are studied by Animasaun et al. [46]. According to the results, the HTG shrinks as the radiation parameter increases. Azam [47] studied the flow of thermally radiative MHD Maxwell nanofluid and CCHF over a surface, considering the effects of a heat source and sink. The rate of heat transfer decreases with increasing values of the radiation parameter.

1.2 Novelty of the study

The key novelties of the proposed work are as follows:

Development of a two-stage exponential time integration scheme tailored explicitly for non-linear PDEs modelling Darcy–Forchheimer–Williamson flow with temperature-dependent thermal conductivity. This formulation has not been addressed in previous literature.
Integrating a sixth-order compact spatial discretization with the proposed time-stepping method enables enhanced resolution and accuracy with fewer grid points.
Fully explicit structure that avoids linearization and iterative solvers, offering both computational efficiency and ease of implementation, particularly valuable for complex boundary layer and porous media problems.
Stability and convergence analysis via the Fourier/von Neumann method, providing theoretical assurance alongside empirical validation.
A comprehensive parametric study was conducted using the proposed high-accuracy framework to demonstrate the significant effects of Weissenberg number, Forchheimer number, and variable thermal conductivity on the fluid behaviour.

Efficient, precise, and stable numerical techniques are essential for solving the governing partial differential equations that describe such flows. Classical numerical methods fail to balance computational cost and accuracy, especially for non-linear, time-dependent partial differential equations. Stability and error reduction require a powerful numerical method to discretize time- and space-dependent variables. This study proposes a new computational method for solving the Darcy–Forchheimer flow equations of Williamson fluids with temperature-dependent thermal conductivity. A unique time-stepping approach provides second-order time precision, while a sixth-order compact finite difference scheme handles spatial discretization. Under specific conditions, the Fourier series or von Neumann stability analysis shows the conditional stability of the system. A popular time-stepping method, the second-order Runge–Kutta method, is compared to the proposed method. Simulations of complex flow phenomena show that the new technique is more accurate and trustworthy due to decreased errors. The study also modifies the planned strategy to improve performance.

The suggested computing scheme is described in Section 2, which also details the spatial and temporal discretization techniques. An examination of the scheme’s robustness is provided in Section 3. The problem’s dimensionless formulation and the governing equations are laid out in Section 4. Section 5 highlights the benefits of the suggested methodology and gives numerical results, including comparisons with the Runge–Kutta method. The results and recommendations for further study are presented in Section 6.

2 Numerical methodology

The numerical scheme will be developed using two time levels. The initial stage of the scheme is explicit and determines the solution at an arbitrary time level. In contrast, the subsequent stage serves as a corrector, identifying the solution at the next time level. To construct the scheme, the whole domain is divided into finite subintervals, and the solution will be found at each endpoint of every subinterval. Let Δ x be the spatial step size, and Δ t be the temporal step size. The differential equation for this case can be written as

(1) ∂ h ∂ t = G h , ∂ h ∂ y , ∂ 2 h ∂ y 2 .

This represents a general parabolic partial differential equation modelling time-dependent fluid behaviour, where h could represent velocity, temperature, or concentration, and G ( ∙ ) includes advection, diffusion, and non-linear source terms, common in fluid and thermal transport problems. This form arises in modelling non-Newtonian fluid flow (e.g. polymer or blood) through porous media, heat transfer where conductivity changes with temperature, and reactive transport in environmental or industrial porous systems.

Under the specified initial and boundary conditions,

(2) h ( 0 , y ) = h 1 , h ( t , 0 ) = h 2 ( t ) , h ( t , L ) = h 3 ( t ) ,

where h 1 is constant and h 2 and h 3 are functions of t . These are needed to start the simulation and define physical constraints, such as prescribed inlet conditions, wall temperatures, or fluid velocities at boundaries, and time-dependent conditions to simulate oscillating environments (e.g. pulsatile flow, thermal cycling).

2.1 First stage – modified exponential integrator

The first stage is a modified form of the exponential time integrator scheme:

(3) h ̅ i n + 1 = h i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 ∂ h ∂ t i n + 0.05 h i n ,

where h ̅ i n + 1 finds the solution at an arbitrary time level and spatial grid point i . This stage uses an explicit exponential time integrator to predict the solution at an intermediate point. It efficiently handles stiff terms in fluid equations (e.g. strong gradients near boundary layers) and avoids the high computational cost of implicit schemes.

2.2 Second stage – corrector via explicit Runge–Kutta form

The second stage of the scheme is the corrector and finding the solution at the next time level n + 1 . The corrector stage is the explicit Runge–Kutta scheme of the form

(4) h i n + 1 = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) ∂ h ̅ ∂ t i n + 1 .

It corrects the predicted solution and advances it in time t n + 1 and uses explicit evaluation; no non-linear solvers are needed. It enables second-order accuracy in time.

Subsequently, examine the Taylor series expansion for h i n + 1 as follows:

(5) h i n + 1 = h i n + Δ t ∂ h ∂ t i n + ( Δ t ) 2 2 ∂ 2 h ∂ t 2 i n + O ( ( Δ t ) 3 ) .

It ensures accuracy and predictable convergence behaviour and matches numerical and analytical solution orders.

Inserting the Taylor series expansion (5) into Eq. (4) yields

(6) h i n + Δ t ∂ h ∂ t i n + ( Δ t ) 2 2 ∂ 2 h ∂ t 2 i n = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) ∂ h ̅ ∂ t i n + 1 .

Since Eq. (6) contains terms that require solution at arbitrary time levels, using the first stage (3) in Eq. (6) yields

(7) h i n + Δ t ∂ h ∂ t i n + ( Δ t ) 2 2 ∂ 2 h ∂ t 2 i n = a h i n + b h i n e 0.05 Δ t + b ( 1 − e − 0.05 Δ t ) 0.05 ∂ h ∂ t i n + 0.05 h i n + c ( e Δ t − 1 ) e − 0.05 Δ t ∂ h ∂ t i n + ( 1 − e − 0.05 Δ t ) 0.05 ∂ 2 h ∂ t 2 i n + 0.05 ∂ h ∂ t i n .

Now, comparing the coefficients of h i n , ∂ h ∂ t i n and ∂ 2 h ∂ t 2 i n on both sides of Eq. (7) yields

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

The proposed scheme for Eq. (10) can be expressed as

(9) h ̅ i n + 1 = h i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 { G i n + 0.05 h i n }

and

(10) h i n + 1 = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) G ̅ i n + 1 ,

where G i n = G h i n , ∂ h ∂ y i n , ∂ 2 h ∂ y 2 i n and G ̅ i n + 1 = G h ̅ i n + 1 , ∂ h ̅ ∂ y i n + 1 , ∂ 2 h ̅ ∂ y 2 i n + 1 .

Let G = α 1 h i n + α 2 ∂ h ∂ y i n + α 3 ∂ 2 h ∂ y 2 i n in Eq. (1), where α 1 is a reaction or damping constant, α 2 is the advection, like transport, and α 3 is the diffusion term. Then, the proposed scheme takes the form:

(11) h ̅ i n + 1 = h i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 α 1 h i n + α 2 ∂ h ∂ y i n + α 3 ∂ 2 h ∂ y 2 i n + 0.05 h i n ,

(12) h i n + 1 = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) α 1 h ̅ i n + 1 + α 2 ∂ h ̅ ∂ y i n + 1 + α 3 ∂ 2 h ̅ ∂ y 2 i n + 1 .

2.3 Spatial discretization – Compact schemes

The compact scheme proposes space discretization. Therefore, by applying compact spatial discretization, Eqs (11) and (12) can be written as

(13) h ̅ i n + 1 = h i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 { α 1 h i n + α 2 P 1 − 1 Q 1 h i n + α 3 P 2 − 1 Q 2 h i n + 0.05 h i n } ,

(14) h i n + 1 = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) { α 1 h ̅ i n + 1 + α 2 P 1 − 1 Q 1 h ̅ i n + 1 + α 3 P 2 − 1 Q 2 h ̅ i n + 1 } .

In Eqs (13) and (14), P 1 , Q 1 , P 2 , and Q 2 are matrices constructed from the coefficients of the following equations:

(15) β 1 h ′ | i − 1 n + h ′ | i n + β 1 h ′ | i + 1 n = c 0 ( h i + 1 n − h i − 1 n ) 2 Δ y + c 1 ( h i + 2 n − h i − 2 n ) 4 Δ y ,

(16) β 1 h ″ | i − 1 n + h ″ | i n + β 2 h ″ | i + 1 n = c 2 ( h i + 1 n − 2 h i n + h i − 1 n ) ( Δ y ) 2 + c 3 ( h i + 2 n − 2 h i n + h i − 2 n ) 4 ( Δ y ) 2 ,

where c 0 = 2 3 ( 2 + β 1 ) , c 1 = 1 3 ( 4 β 1 − 1 ) , c 2 = 4 3 ( 1 − β 2 ) , c 3 = 1 3 ( 10 β 2 − 1 ) .

Using sixth-order compact finite difference schemes gives high resolution of gradients near boundary layers (common in heat/mass transfer), reduces numerical dispersion, and improves stability.

3 Stability analysis

One of the established stability criteria for assessing the stability conditions of finite difference schemes is the Von Neumann stability analysis, sometimes referred to as the Fourier series stability analysis. The criterion is based on transforming the given difference equation into trigonometric equations that further produce stability conditions. First, the given differential equation is discretized by the finite difference scheme, and transformations are employed. For the present case, the transformations can be expressed as

(17) P 1 e iI ψ = β 1 e ( i + 1 ) I ψ + e iI ψ + β 1 e ( i − 1 ) I ψ ,

(18) Q 1 e iI ψ = c 0 ( e ( i + 1 ) I ψ − e ( i − 1 ) I ψ ) 2 Δ y + c 1 ( e ( i + 2 ) I ψ − e ( i − 2 ) I ψ ) 4 Δ y ,

(19) P 2 e iI ψ = β 2 e ( i + 1 ) I ψ + e iI ψ + β 2 e ( i − 1 ) I ψ ,

(20) Q 2 e iI ψ = c 2 ( e ( i + 1 ) I ψ − 2 e iI ψ + e ( i − 1 ) I ψ ) ( Δ y ) 2 + c 3 ( e ( i + 2 ) I ψ − 2 e iI ψ + e ( i − 2 ) I ψ ) 4 ( Δ y ) 2 ,

where I = − 1 is an imaginary number.

The first stage of the proposed scheme, after employing transformations (17)–(20), gives

(21) h ̅ i n + 1 = h i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 α 1 h i n + α 2 ( 4 c 0 I sin ψ + c 1 I sin 2 ψ ) 2 Δ y ( 2 β 1 cos ψ + 1 ) h i n + α 3 8 c 2 ( cos ψ − 1 ) + 2 c 3 ( cos 2 ψ − 1 ) 4 ( Δ y ) 2 ( 2 β 2 cos ψ + 1 ) h i n + 0.05 h i n .

Rewrite Eq. (21) as

(22) h ̅ i n + 1 = ( δ 1 + i δ 2 ) h i n ,

where δ 1 = e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 α 1 + α 3 8 c 2 ( cos ψ − 1 ) + 2 c 3 ( cos 2 ψ − 1 ) 4 ( Δ y ) 2 ( 2 β 2 cos ψ + 1 ) + 0.05 ,

δ 2 = ( 1 − e − 0.05 Δ t ) 0.05 α 2 ( 4 c 0 I sin ψ + c 1 I sin 2 ψ ) 2 Δ y ( 2 β 1 cos ψ + 1 ) .

By employing transformations (17)–(20) into second stage of the scheme, it yields

(23) h i n + 1 = a h i n + b h ̅ i n + 1 + c ( e Δ t − 1 ) α 1 h ̅ i n + 1 + α 2 ( 4 c 0 I sin ψ + c 1 I sin 2 ψ ) 2 Δ y ( 2 β 1 cos ψ + 1 ) h ̅ i n + 1 + α 3 8 c 2 ( cos ψ − 1 ) + 2 c 3 ( cos 2 ψ − 1 ) 4 ( Δ y ) 2 ( 2 β 2 cos ψ + 1 ) h ̅ i n + 1 .

Eq. (23) can be reformulated as follows:

(24) h i n + 1 = δ 3 h i n + ( δ 4 + I δ 5 ) h ̅ i n + 1 ,

where δ 3 = a , δ 4 = b + c ( e Δ t − 1 ) α 1 + α 3 8 c 2 ( cos ψ − 1 ) + 2 c 3 ( cos 2 ψ − 1 ) 4 ( Δ y ) 2 ( 2 β 2 cos ψ + 1 ) and δ 5 = c ( e Δ t − 1 ) α 2 ( 4 c 0 I sin ψ + c 1 I sin 2 ψ ) 2 Δ y ( 2 β 1 cos ψ + 1 ) .

Inserting Eq. (22) into Eq. (24),

(25) h i n + 1 = ( δ 6 + I δ 7 ) h i n ,

where δ 6 = δ 3 + δ 1 δ 4 − δ 5 δ 2 and δ 7 = δ 4 δ 2 + δ 5 δ 1 .

The amplification factor for this case can be written as

(26) h i n + 1 h i n 2 = δ 6 2 + δ 7 2 ≤ 1 .

If the proposed scheme satisfies inequality (26), it will remain stable; otherwise, the solution will diverge. Therefore, for a converged solution, the step size and parameters involved in the given differential equation must satisfy inequality (26).

The stability criterion for the scalar partial differential equation was provided. The convergence of the partial differential equations will now be presented. To do so, consider the vector–matrix equation of the convection–diffusion equation as

(27) ∂ g ∂ t = C 1 g + C 2 ∂ g ∂ y + C 3 ∂ 2 g ∂ y 2 ,

where g is a vector and C 1 , C 2 , and C 3 are matrices.

Applying the proposed scheme to Eq. (27) for time discretization and employing the compact scheme to discretize space-dependent terms in Eq. (27),

(28) g ̅ i n + 1 = g i n e − 0.05 Δ t + 1 − e − 0.05 Δ t 0.05 { C 1 g i n + C 2 P 1 − 1 Q 1 g i n + C 3 P 2 − 1 Q 2 g i n + 0.05 g i n } ,

(29) g i n + 1 = a g i n + b g ̅ i n + 1 + c ( e Δ t − 1 ) { C 1 g ̅ i n + 1 + C 2 P 1 − 1 Q 1 g ̅ i n + 1 + C 3 P 2 − 1 Q 2 g ̅ i n + 1 } .

Theorem

The proposed computational scheme in a time and compact scheme in space for Eq. (27) converges.

Proof

In the beginning, for the proof of this theorem, consider the exact scheme for Eq. (27) as

(30) G ̅ i n + 1 = G i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 { C 1 G i n + C 2 P 1 − 1 Q 1 G i n + C 3 P 2 − 1 Q 2 G i n + 0.05 G i n } ,

(31)□ G i n + 1 = a G i n + b G ̅ i n + 1 + c ( e Δ t − 1 ) { C 1 G ̅ i n + 1 + C 2 P 1 − 1 Q 1 G ̅ i n + 1 + C 3 P 2 − 1 Q 2 G ̅ i n + 1 } .

Now, subtracting the first stage of the proposed scheme (28) from the first stage of the exact scheme that yields G i n − g i n = E i n , etc.,

(32) E ̅ i n + 1 = E i n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 { C 1 E i n + C 2 P 1 − 1 Q 1 E i n + C 3 P 2 − 1 Q 2 E i n + 0.05 E i n } .

By taking the norm ‖ ∙ ‖ ∞ on both sides of Eq. (32),

(33) E ̅ n + 1 = E n e − 0.05 Δ t + ( 1 − e − 0.05 Δ t ) 0.05 { ‖ C 1 ‖ ∞ E n + ‖ C 2 P 1 − 1 Q 1 ‖ ∞ E n + ‖ C 3 P 2 − 1 Q 2 ‖ ∞ E n + 0.05 E n } .

Inequality (33) can be rewritten as

(34) E ̅ n + 1 ≤ λ 1 E n ,

where λ 1 = e − 0.05 Δ t ( 1 − e − 0.05 Δ t ) 0.05 { ‖ C 1 ‖ ∞ E n + ‖ C 2 P 1 − 1 Q 1 ‖ ∞ E n + ‖ C 3 P 2 − 1 Q 2 ‖ ∞ E n + 0.05 E n } .

The current result is derived by deducting the second stage of the suggested scheme (29) from the second stage of the exact scheme (31):

(35) E i n + 1 = a E i n + b E ̅ i n + 1 + c ( e Δ t − 1 ) { C 1 E ̅ i n + 1 + C 2 P 1 − 1 Q 1 E ̅ i n + 1 + C 3 P 2 − 1 Q 2 E ̅ i n + 1 } .

By taking the norm ‖ ∙ ‖ ∞ on both sides of Eq. (35),

(36) E n + 1 = | a | E n + | b | E ̅ n + 1 + | c | | e Δ t − 1 | { ‖ C 1 ‖ ∞ + ‖ C 2 P 1 − 1 Q 1 ‖ ∞ + ‖ C 3 P 2 − 1 Q 2 ‖ ∞ } E ̅ n + 1 .

By rewriting inequality (36) as

(37) E n + 1 ≤ λ 2 E n + λ 3 E ̅ n + 1 ,

where μ 3 = b + | c | | e Δ t − 1 | { ‖ C 1 ‖ ∞ + ‖ C 2 P 1 − 1 Q 1 ‖ ∞ + ‖ C 3 P 2 − 1 Q 2 ‖ ∞ } and λ 2 = | a | .

Using inequality (34) in inequality (37) yields

E n + 1 ≤ λ 2 E n + λ 3 E ̅ n + 1 = E n + 1 ≤ λ 2 E n + λ 3 λ 1 E n ,

(38) E n + 1 ≤ ( λ 2 + λ 1 λ 3 ) E n .

Now, rewrite inequality (38) as

(39) E n + 1 ≤ λ 4 E n + R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) .

Let n = 0 in inequality (39). It yields

(40) E 1 ≤ λ 4 E 0 + R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) .

Due to the exact initial condition E 0 = 0 , inequality (40) yields

(41) E 1 ≤ R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) .

Let n = 1 in inequality (39). It gives

(42) E 2 ≤ λ 4 E 1 R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) ≤ ( 1 + λ 4 ) R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) .

If this is continued, then for finite n ,

(43) E n ≤ ( 1 + λ 4 + … + λ 4 n − 1 ) R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) = 1 − λ 4 n 1 − λ 4 R ( O ( ( Δ t ) 2 , ( Δ y ) 6 ) ) .

By applying limit when n → ∞ , the series 1 + λ 4 + … + λ 4 n + … becomes an infinite geometric series that converges when | λ 4 | < 1 .

4 Mathematical formulation

Consider a modelling of unsteady, laminar flow of a non-Newtonian (Williamson) fluid near a vertical flat plate embedded in a porous medium. This setup models cooling systems, polymer extrusion, oil recovery, and biofluid transport, where non-Newtonian properties, temperature variation, and porosity significantly influence flow behaviour.

Geometry: Let x * be the vertical direction along the plate, y * be normal to the plate (direction of flow development), and driving forces be the buoyancy due to temperature and concentration gradients (i.e. thermal and solutal convection). The temperature and concentration gradient drive the flow in the fluid. Let the temperature and concentration away from the plate be less than those at the plate. At t * = 0 , fluid is at ambient conditions T ∞ , C ∞ , while for t * > 0 , oscillatory heating begins. The considered model also uses the effect of Darcy–Forchheimer’s characteristics. Figure 1 depicts the physical arrangement of the problem with unstable, mixed convective flow of a Williamson fluid over a vertical flat plate situated within a porous material.

Figure 1

Geometry of the problem.

Under this scenario, the governing equation of this model can be expressed as

(44) ∂ u * ∂ t * = ν ∂ 2 u * ∂ y * 2 + 2 ν Γ ∂ u * ∂ y * ∂ 2 u * ∂ y * 2 + g β T ( T − T ∞ ) + g β c ( C − C ∞ ) − σ B o 2 ρ u * − ν k p u * − b ̅ u * 2 ,

(45) ∂ T ∂ t * = 1 ρ c p ∂ ∂ y * k ∞ ( 1 + ϵ 1 θ ) ∂ T ∂ y * + ν c p k p μ * 2 + b c p ̅ μ * 3 + q ″ ′

(46) ∂ C ∂ t * = D ∂ 2 C ∂ y * 2 − k 1 ( C − C ∞ ) .

Here, we explain these governing equations in terms of their physical interpretation using a world scenario.

Momentum (Eq. ( 44 )): ∂ u * ∂ t * represents the time derivative of velocity; unsteady flow also describes how the horizontal velocity u * evolves over time. ν ∂ 2 u * ∂ y * 2 represents viscous diffusion. The standard Newtonian viscous term describes how momentum diffuses from the wall into the fluid (boundary layer development). 2 ν Γ ∂ u * ∂ y * ∂ 2 u * ∂ y * 2 represents the Williamson fluid elasticity term, where Γ stands for fluid relaxation time; this non-linear term accounts for shear thinning and viscoelastic effects, and also reflects the non-Newtonian behaviour of Williamson fluids. g β T ( T − T ∞ ) accounts for thermal buoyancy force, as T > T ∞ the heated fluid becomes less dense and increases, and drives natural convection along the vertical plate. g β c ( C − C ∞ ) represents Solutal (concentration) buoyancy force, analogous to thermal buoyancy but due to concentration differences it is also significant in double-diffusive convection. − σ B o 2 ρ u * accounts for magnetic (Lorentz) force, and acts as a resisting force in the presence of a magnetic field B o . − ν k p u * represents Darcy’s resistance and models linear drag due to fluid moving through porous media, where k p is the permeability of the medium. − b ̅ u * 2 accounts for Forchheimer resistance and models inertial (non-linear) drag in porous media, important at higher velocities.

Energy (Eq. ( 45 )): ∂ T ∂ t * represents the temperature change rate with time, unsteady (transient) heating/cooling. 1 ρ c p ∂ ∂ y * k ∞ ( 1 + ϵ 1 θ ) ∂ T ∂ y * represents the conduction term, which includes temperature-dependent thermal conductivity, where ϵ 1 is the sensitivity parameter when ϵ 1 = 0 , and conductivity is constant. ν c p k p μ * 2 + b c p ̅ μ * 3 + q ″ ′ accounts for dissipation terms, representing viscous and Forchheimer-related heating; it is important in porous flows, where q ″ ′ is a space and temperature-dependent internal heat generation and is defined as q ″ ′ = k ∞ u R ρ x * ν c p ( A * u ( T w − T ∞ ) + B * ( T − T ∞ ) ) , where A * denotes the space-dependent heat generation coefficient, which accounts for the influence of velocity and spatial variation on internal heating near the boundary, and B * denotes the temperature-dependent heat generation coefficient, which modulates the heat generation based on the deviation of local temperature from the ambient state. These parameters are commonly used in modelling radiative and reactive heat generation processes within porous or non-homogeneous media.

Mass diffusion (Eq. ( 46 )): Eq. (46) describes diffusion with chemical reaction, where D is the mass diffusivity, k 1 is the reaction rate parameter, C ∞ is the ambient concentration, D ∂ 2 C ∂ y * 2 represents the Fickian diffusion, and − k 1 ( C − C ∞ ) accounts for chemical consumption or decay subject to initial and boundary conditions,

(47) u * = 0 , T = T ∞ , C = C ∞ when t * = 0 u * = 0 , T = T ∞ + ϵ 2 ( T w − T ∞ ) cos w * t * , C = C ∞ + ϵ 2 ( C w − C ∞ ) cos w * t * when y * = 0 u * → 0 , T → 0 , C → 0 when y * → ∞ .

Initial and boundary conditions (Eq. (47)): Initial state t * = 0 system is at rest with uniform temperature and concentration. At the wall, y * = 0 , the velocity is zero (no-slip condition), and temperature and concentration fluctuate periodically in time (cosine wave), simulating oscillating surface conditions. Far from the wall, y * → ∞ , field variables return to the ambient state.

Let B o be the strength of the magnetic field applied transversely to the plate, ν be the kinematic viscosity, σ be the electrical conductivity, β T and β c be, respectively, the coefficients of thermal and solutal convection, k 1 be the reaction rate, b ̅ be the non-Darcian parameter, g be the gravity, ρ be the density of the fluid, and k p denotes the permeability constant.

Dimensionless transformation (Eq. ( 48 )): Now, using the transformations

(48) y = y * L R , u = u * u R , w = t R w * , t = t * t R , θ = T − T ∞ T w − T ∞ , ϕ = C − C ∞ C w − C ∞ ,

where u R = ( ν g β T Δ T ) 1 3 is the reference velocity due to natural convection, L R = g β T Δ T ν 2 − 1 3 is the boundary layer thickness, and t R = ( g β T Δ T ) − 2 3 ν 1 3 is the time scale based on thermal buoyancy. This will simplify the equations and reduce the number of parameters.

Dimensionless PDE system: The governing Eqs (44)–(47) are reduced to

(49) ∂ u ∂ t = ∂ 2 u ∂ y 2 + W e ∂ 2 u ∂ y 2 ∂ u ∂ y − M + 1 D a u − F s u 2 + θ + NC ,

(50) ∂ θ ∂ t = 1 P r ( 1 + ϵ 1 θ ) ∂ 2 θ ∂ y 2 + ϵ 1 ∂ θ ∂ y 2 + ϵ P r ( A * u + B * θ ) + E c D a u 2 + E c F s u 3 ,

(51) ∂ ϕ ∂ t = 1 S c ∂ 2 ϕ ∂ y 2 − k c ϕ ,

subject to the dimensionless initial and boundary conditions

(52) u = 0 , θ = 0 , ϕ = 0 for t = 0 u = 0 , θ = ϵ 2 cos wt , ϕ = ϵ 2 cos wt for y = 0 u → 0 , θ → 0 , ϕ → 0 for y → ∞ ,

where N represents the buoyancy ratio, D a is Darcy’s number, M is a magnetic parameter, E c is Eckert’s number, P r is the Prandtl number, W e is the Weissenberg number, S c is the Schmidt number, ϵ 2 is the dimensionless parameter, F S represents the Forchheimer number, k c is a dimensionless reaction rate parameter, and these parameters are given as

N = β c ( C w − C ∞ ) β T ( T w − T ∞ ) , D a = k p ν t R , M = t R σ B o 2 ρ , E c = u R 2 c p ( T w − T ∞ ) , P r = ν k ∞ ρ c p , W e = Γ u R 2 L R , ϵ 2 = t R u R x * , F s = b ̅ u R t R , k c = t R k c , S c = ν D .

The local Nusselt number N uL and local Sherwood number S h L are defined as

(53) N u L = q w L R k ∞ ( T w − T ∞ ) , S h L = j w L R D ( C w − C ∞ ) ,

where q w = − k ∞ ∂ T ∂ y * y * = 0 is the wall heat flux and j w = − D ∂ C ∂ y * y * = 0 is the wall mass flux.

Using Eq. (48), the dimensionless local Nusselt and Sherwood numbers are given as

(54) N u L = − ∂ θ ∂ y y = 0 S h L = − ∂ ϕ ∂ y y = 0 .

5 Results and discussions

We conducted a comprehensive simulation research with the subsequent aims.

5.1 Demonstration of the proposed numerical scheme

As developed in Section 5, the proposed numerical scheme is applied to solve the governing equations for the Darcy–Forchheimer flow of Williamson fluid with temperature-dependent thermal conductivity. The scheme preserves the accuracy and stability of the solution by employing a second-order explicit time-stepping technique combined with a sixth-order compact scheme for spatial discretization. The explicit nature of the scheme eliminates the need for linearization or iterative solvers, significantly reducing the computational time. Numerical results confirm the scheme’s conditional stability, validated using Fourier series or von Neumann stability analysis. The positivity of the solution is also maintained, ensuring the physical relevance of the results.

5.2 Stability and convergence of the scheme

The stability of the proposed scheme is demonstrated under specific conditions on step sizes and parameters. The scheme is consistent, verified by Taylor series expansion, and convergent for linear partial differential equations, as guaranteed by the Lax equivalence theorem. Numerical experiments show that the scheme performs reliably for a wide range of parameter values, accurately capturing the flow’s dynamics. This is further confirmed by the attainment of steady-state solutions under stable conditions.

5.3 Velocity profile analysis

Figure 2 illustrates the variation of the velocity profile for different values of the Weissenberg number W e while keeping the other parameters constant. The Weissenberg number is a dimensionless quantity that characterizes elastic effects in non-Newtonian (viscoelastic) fluids. It is defined as W e = relaxation time × characteristic velocity characteristic length and quantifies the ratio of elastic to viscous forces. As W e increases from 0.1 to 0.9, the velocity profile becomes slightly higher near the wall and gradually decays into the free stream region. The solid line for W e = 0.1 shows the steepest decline, while the dotted line for W e = 0.9 indicates a more extended boundary layer. A higher W e implies a stronger elastic memory effect of the fluid particles. The result is a thicker momentum boundary layer, which is consistent with the viscoelastic behaviour of Williamson fluids.

$Figure 2 Variation of the Weissenberg number on the velocity profile using F s = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , A * = 0.1 , {F}_{{\rm{s}}}=0.1,{D}_{{\rm{a}}}=7,N=0.1,M=0.1,{A}^{* }=0.1, B * = 0.1 , P r = 0.9 , {B}^{* }=0.1,{P}_{{\rm{r}}}=0.9, γ = 0.1 , S c = 0.9 , \gamma =0.1,{S}_{{\rm{c}}}=0.9, ε = 0.1 , and ε 1 = 0.1 \varepsilon =0.1,{\rm{and}}{\varepsilon }_{1}=0.1 .$

Figure 2

Variation of the Weissenberg number on the velocity profile using F s = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , γ = 0.1 , S c = 0.9 , ε = 0.1 , and ε 1 = 0.1 .

Figure 3 depicts the variation of the velocity profile for different values of the Forchheimer number F S , while keeping all other parameters constant. The velocity profile decreases as the Forchheimer number F S increases from 0.1 to 0.9. The boundary layer becomes thinner, meaning the fluid loses momentum closer to the wall. The solid line F s = 0.1 shows the highest velocity, dashed line F s = 0.5 shows a moderate reduction, and the dotted line F s = 0.9 shows the strongest damping in velocity. With increasing F s , the fluid experiences greater inertial resistance, especially at higher velocities near the wall. This enhanced resistance suppresses the flow, causing a sharper velocity decay away from the wall. The impact is especially evident in porous media where non-linear drag effects dominate at higher flow intensities or permeability gradients.

$Figure 3 Variation of the Forchheimer number on the velocity profile using W e = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=7,N=0.1,M=0.1,{A}^{* }=0.1,{B}^{* }=0.1,{P}_{{\rm{r}}}=0.9, γ = 0.1 , S c = 0.9 , \gamma =0.1,{S}_{{\rm{c}}}=0.9, ε = 0.1 , and ε 1 = 0.1 \varepsilon =0.1,{\rm{and}}{\varepsilon }_{1}=0.1 .$

Figure 3

Variation of the Forchheimer number on the velocity profile using W e = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , γ = 0.1 , S c = 0.9 , ε = 0.1 , and ε 1 = 0.1 .

Figure 4 demonstrates the effect of the buoyancy parameter N on the velocity profile, with other parameters held constant. The buoyancy parameter N (also called the buoyancy ratio) is defined as N = β c ( C w − C ∞ ) β T ( T w − T ∞ ) . It represents the solutal (concentration-driven) buoyancy’s relative effect on thermal (temperature-driven) buoyancy in mixed convective flow. The velocity profile increases with the higher buoyancy parameter N ; solid lines: N = 0.1 , dashed line: N = 0.5 , and dotted line: N = 0.9 . Higher values of N lead to slower velocity decay, indicating a thicker momentum boundary layer. As N increases, the influence of solutal buoyancy becomes more dominant relative to thermal buoyancy. This additional upward driving force (from concentration gradients) boosts the fluid velocity. As a result, the momentum in the fluid increases, especially near the plate. The fluid retains velocity farther from the wall, leading to an extended velocity profile.

$Figure 4 Variation of the buoyancy parameter on the velocity profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5,{F}_{{\rm{s}}}=0.1,M=0.1, A * = 0.1 , {A}^{* }=0.1, B * = 0.1 , P r = 0.9 , γ = 0.1 , {B}^{* }=0.1,{P}_{{\rm{r}}}=0.9,\gamma =0.1, S c = 0.9 , ε = 0.1 , and ε 1 = 0.1 {S}_{{\rm{c}}}=0.9,\varepsilon =0.1,{{and}\varepsilon }_{1}=0.1 .$

Figure 4

Variation of the buoyancy parameter on the velocity profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , γ = 0.1 , S c = 0.9 , ε = 0.1 , and ε 1 = 0.1 .

5.4 Impact of thermal conductivity and other parameters on temperature and concentration profiles

Effect of thermal conductivity parameter on the temperature profile: Figure 5 depicts the variation of the temperature profile for different values of the small parameter of thermal conductivity ϵ 1 while keeping other parameters constant. With an increase in the tiny parameter ϵ 1 , the temperature profile increases. After the temperature profile peaks, a steady-state value is attained farther from the boundary. As the value of ϵ 1 increases, the temperature profile becomes flatter, suggesting better heat distribution. The greater the value of ϵ 1 , the fluid’s thermal conductivity increases more dramatically as temperature increases. As ϵ 1 increases, the heat transfer within the fluid increases, and the heat transfer becomes more effective, resulting in an overall increase in the temperature profile. The enhanced thermal conductivity leads to better heat diffusion, which explains the elevated and flattened temperature profiles for larger ϵ 1 .

$Figure 5 Variation of the small parameter of thermal conductivity on the temperature profile using W e = 0.1 , D a = 5 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5, F s = 0.1 , M = 0.1 , A * = 0.1 , {F}_{{\rm{s}}}=0.1,M=0.1,{A}^{* }=0.1, B * = 0.1 , P r = 0.9 , γ = 0.1 , {B}^{* }=0.1,{P}_{{\rm{r}}}=0.9,\gamma =0.1, S c = 0.9 , ε = 0.1 , and N = 0.1 {S}_{{\rm{c}}}=0.9,\varepsilon =0.1,{\rm{and}}N=0.1 .$

Figure 5

Variation of the small parameter of thermal conductivity on the temperature profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , γ = 0.1 , S c = 0.9 , ε = 0.1 , and N = 0.1 .

Effect of reaction rate parameter on the concentration profile: Figure 6 illustrates the reaction rate parameter γ on the concentration profile, with other constant parameters. The concentration profile decreases with increasing values of the reaction rate parameter γ . Higher values of γ cause a noticeable reduction in the peak concentration, as well as a quicker drop-off in the concentration further from the boundary. The concentration profile stabilizes to a steady-state value faster for higher γ . The reaction rate parameter γ reflects the intensity of chemical reactions within the fluid. A higher γ indicates faster reaction rates. As γ increases, reactants are consumed more rapidly during the chemical reactions, leading to a lower concentration of reactants in the fluid. This explains the overall decline in the concentration profile for larger γ . The reduced concentration near the boundary is due to the faster depletion of reactants in the reaction zone. At the same time, the stabilization further away represents the diffusion-dominated region where reaction effects are less pronounced.

$Figure 6 Variation of the reaction rate parameter on the concentration profile using W e = 0.1 , D a = 5 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5, F s = 0.1 , M = 0.1 , {F}_{{\rm{s}}}=0.1,M=0.1, A * = 0.1 , B * = 0.1 , {A}^{* }=0.1,{B}^{* }=0.1, P r = 0.9 , ε 1 = 0.1 , {P}_{{\rm{r}}}=0.9,{\varepsilon }_{1}=0.1, S c = 0.9 , ε = 0.1 , N = 0.1 {S}_{{\rm{c}}}=0.9,\varepsilon =0.1,N=0.1 .$

Figure 6

Variation of the reaction rate parameter on the concentration profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , P r = 0.9 , ε 1 = 0.1 , S c = 0.9 , ε = 0.1 , N = 0.1 .

5.5 Effect on local Nusselt and Sherwood numbers

Figure 7 shows the variation of the local Nusselt number N u L for different values of the small parameter of thermal conductivity ϵ 1 and Prandtl number P r , while keeping the other parameters constant. The local Nusselt number decreases with an increase in the small parameter of thermal conductivity ϵ 1 . A similar declining trend is observed for increasing values of the Prandtl number P r . Higher ϵ 1 and P r values lead to a more rapid reduction in the local Nusselt number. The local Nusselt number N u L represents the heat transfer rate at the surface. A higher ϵ indicates an increased dependency of thermal conductivity on temperature, which enhances heat diffusion within the fluid, reducing the surface heat transfer rate. The Prandtl number P r characterizes the ratio of momentum diffusivity to thermal diffusivity. Higher P r values imply slower thermal diffusion compared to momentum diffusion, leading to a reduced surface heat transfer rate.

$Figure 7 Variation of the Prandtl number and small parameter of thermal conductivity on the local Nusselt number using W e = 0.1 , D a = 5 , {W}_{{\rm{e}}}=0.1,\hspace{0.25em}{D}_{{\rm{a}}}=5, F s = 0.1 , M = 0.1 , {F}_{{\rm{s}}}=0.1,{M}=0.1, A * = 0.1 , B * = 0.1 , {A}^{* }=0.1,\hspace{0.25em}{B}^{* }=0.1, γ = 0.1 , S c = 0.9 , \gamma =0.1,\hspace{0.25em}{S}_{{\rm{c}}}=0.9, ε = 0.1 , \varepsilon =0.1, and N = 0.1 {\rm{and}}{N}=0.1 .$

Figure 7

Variation of the Prandtl number and small parameter of thermal conductivity on the local Nusselt number using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , γ = 0.1 , S c = 0.9 , ε = 0.1 , and N = 0.1 .

Figure 8 illustrates the variation of the local Sherwood number S h L for different values of the reaction rate parameter γ and Schmidt number S c , while other parameters are kept constant. The local Sherwood number decreases with an increase in the reaction rate parameter γ . Similarly, the local Sherwood number decreases as the Schmidt number S c increases. The reduction is more pronounced for higher values of both γ and S c . The local Sherwood number S h indicates the mass transfer rate at the surface. A higher reaction rate parameter γ means faster chemical reactions, consuming more reactants near the surface and leading to a reduced mass transfer rate. The Schmidt number S c reflects the ratio of momentum diffusivity to mass diffusivity. Higher S c values imply slower mass diffusion than momentum diffusion, which limits the transport of reactants to the surface, reducing the Sherwood number.

$Figure 8 Variation of the Schmidt number and reaction rate parameter on the local Sherwood number using W e = 0.1 , D a = 5 , F s = 0.1 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5,{F}_{{\rm{s}}}=0.1, M = 0.1 , A * = 0.1 , M=0.1,{A}^{* }=0.1, B * = 0.1 , ε 1 = 0.1 , Pr = 0.9 , ε = 0.1 , and N = 0.1 {B}^{* }=0.1,{\varepsilon }_{1}=0.1,{\rm{Pr}}=0.9,\varepsilon =0.1,{\rm{and}}N=0.1 .$

Figure 8

Variation of the Schmidt number and reaction rate parameter on the local Sherwood number using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , ε 1 = 0.1 , Pr = 0.9 , ε = 0.1 , and N = 0.1 .

5.6 Contour plots for velocity, temperature, and concentration

Figure 9 presents a contour plot of the velocity profile for a non-Newtonian Williamson fluid under the given parameters. The velocity distribution exhibits periodic patterns near the surface, represented by alternating high- and low-velocity regions (red and blue areas, respectively). The velocity is higher near the boundary (bottom left corner) and decreases gradually as we move away into the bulk fluid (towards the upper region). The colour bar on the right quantifies the velocity magnitudes, with red representing the maximum velocity and blue the minimum. The contour plot reflects the effect of mixed convective flow near the surface, where oscillatory boundary conditions are applied. Both thermal and solutal buoyancy forces influence the velocity near the surface. The periodic pattern near the bottom represents flow instabilities due to the porous medium’s interaction between buoyancy forces and inertial effects. The Forchheimer term enhances this interaction F s = 0.1 , which accounts for inertial resistance.

$Figure 9 Contour plot for the velocity profile using W e = 0.1 , D a = 5 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5, F s = 0.1 , M = 0.1 , {F}_{{\rm{s}}}=0.1,M=0.1, A * = 0.1 , B * = 0.1 , ε 1 = 0.1 , {A}^{* }=0.1,{B}^{* }=0.1,{\varepsilon }_{1}=0.1, P r = 0.9 , ε = 0.1 , {P}_{{\rm{r}}}=0.9,\varepsilon =0.1, N = 0.1 , N=0.1, S c = 0.9 , and γ = 0.9 {S}_{{\rm{c}}}=0.9,{\rm{and}}\gamma =0.9 .$

Figure 9

Contour plot for the velocity profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , ε 1 = 0.1 , P r = 0.9 , ε = 0.1 , N = 0.1 , S c = 0.9 , and γ = 0.9 .

Figure 10 presents a contour plot of the temperature profile for a non-Newtonian Williamson fluid under the specified parameters. The temperature distribution exhibits periodic patterns near the surface, with alternating regions of higher and lower temperature (red and blue areas, respectively). As one advances deeper into the bulk fluid, the temperature drops from its peak close to the border. As indicated by the periodic behaviour, boundary conditions and buoyant forces influence the oscillatory thermal effects. The red zones, regions of high temperature close to the surface, are caused by heat transfer from the barrier into the fluid. The imposed oscillatory boundary conditions and mixed convective influences cause these regions to oscillate. The dominance of thermal diffusion decreases the temperatures with increasing distance from the boundary, leading to a steady decrease in heat transfer farther from the surface. The fluctuation in thermal conductivity, which is affected by the Prandtl number P r = 0.9 , and the tiny parameter ϵ 1 = 0.1 , as well as the buoyancy-driven flow, is reflected in the periodicity of the contours.

$Figure 10 Contour plot for the temperature profile using W e = 0.1 , D a = 5 , Fs = 0.1 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5,{\rm{Fs}}=0.1, M = 0.1 , A * = 0.1 , M=0.1,{A}^{* }=0.1, B * = 0.1 , {B}^{* }=0.1, ε 1 = 0.1 , {\varepsilon }_{1}=0.1, P r = 0.9 , ε = 0.1 , {P}_{{\rm{r}}}=0.9,\varepsilon =0.1, N = 0.1 , S c = 0.9 , and γ = 0.9 N=0.1,{S}_{{\rm{c}}}=0.9,{\rm{and}}\gamma =0.9 .$

Figure 10

Contour plot for the temperature profile using W e = 0.1 , D a = 5 , Fs = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , ε 1 = 0.1 , P r = 0.9 , ε = 0.1 , N = 0.1 , S c = 0.9 , and γ = 0.9 .

Figure 11 presents a contour plot of the concentration profile for a non-Newtonian Williamson fluid under the specified parameters. The concentration profile exhibits periodic patterns near the surface, with regions of high and low concentrations (red and blue areas, respectively). The concentration is highest near the boundary and gradually decreases as it moves farther into the bulk fluid. The oscillatory patterns near the surface become smoother and diffuse as the distance from the boundary increases. The high-concentration zones near the surface (red areas) indicate regions with a higher concentration of solutes. Boundary conditions and reactive effects influence these areas. The gradual decrease in concentration farther from the surface is due to diffusion and consumption of reactants during chemical reactions. The Schmidt number S c = 0.9 plays a key role, indicating the relative importance of mass diffusion versus momentum diffusion. The oscillatory nature of the contours near the boundary reflects the imposed oscillatory boundary conditions and the coupling between mass transfer and flow dynamics.

$Figure 11 Contour plot for the concentration profile using W e = 0.1 , D a = 5 , {W}_{{\rm{e}}}=0.1,{D}_{{\rm{a}}}=5, F s = 0.1 , M = 0.1 , {F}_{{\rm{s}}}=0.1,M=0.1, A * = 0.1 , B * = 0.1 , {A}^{* }=0.1,{B}^{* }=0.1, ε 1 = 0.1 , P r = 0.9 , ε = 0.1 , {\varepsilon }_{1}=0.1,{P}_{{\rm{r}}}=0.9,\varepsilon =0.1, N = 0.1 , S c = 0.9 , and γ = 0.9 N=0.1,{S}_{{\rm{c}}}=0.9,{\rm{and}}\gamma =0.9 .$

Figure 11

Contour plot for the concentration profile using W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , A * = 0.1 , B * = 0.1 , ε 1 = 0.1 , P r = 0.9 , ε = 0.1 , N = 0.1 , S c = 0.9 , and γ = 0.9 .

5.7 Comparison with standard schemes

To evaluate the efficiency of the proposed scheme, it is compared with the second-order Runge–Kutta method. The modified scheme, employing an alternative formulation for the second stage, exhibits lower errors, as shown in Table 1. If the second stage, Eq. (3) of the proposed scheme, is replaced with

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

Table 1

Comparison of the proposed and existing second-order Runge–Kutta method for Stokes’ first problem using N x ( grid points ) = 50 and t f = 1

N t	L 2 error
	Proposed		Second-order Runge–Kutta
	Central	Compact	Central	Compact
250	0.0858	0.0747	0.0876	0.0770
300	0.0754	0.0732	0.0768	0.0751
350	0.0695	0.0730	0.0705	0.0746
400	0.0669	0.0734	0.0677	0.0748
450	0.0667	0.0743	0.0674	0.0755
500	0.0680	0.0753	0.0685	0.0763

the resultant strategy yields fewer mistakes than the second-order Runge–Kutta approach. Table 1 presents a comparison of the L 2 error between the suggested technique and the established second-order Runge–Kutta approach for Stokes’ first issue. The errors are computed using two spatial discretization techniques: central difference and compact scheme. The results are presented for different time-step sizes corresponding to N t (number of time steps), with the grid points in the spatial domain fixed at N x = 50 , and the final time t f = 1 . The proposed scheme consistently produces lower L 2 errors compared to the existing second-order Runge–Kutta method for both central difference and compact schemes, highlighting its superior accuracy.

5.8 Impact of the numerical method on parameter sensitivity

The proposed numerical method enhances the accuracy of capturing the effects of dimensionless parameters by combining a second-order time integrator with a sixth-order compact spatial discretization. This high-resolution approach allows for the precise prediction of how key dimensionless numbers, such as the Weissenberg W e , Forchheimer F s , and Prandtl number P r influence flow, thermal, and concentration fields, without requiring linearization or simplification of non-linear terms. Table 2 provides a structured overview of all the parametric cases studied in this research work. It lists the key dimensionless numbers varied during simulations, such as the Weissenberg number, Forchheimer number, thermal conductivity parameter, their tested ranges, and the associated fixed parameter values. This systematic classification helps isolate the effect of each parameter on the fluid flow behaviour, heat, and mass transfer.

Table 2

Summary of controlling parameter ranges

Case No.	Investigated parameter	Symbol	Range/values tested	Fixed values for other parameters	Related figure
1	Weissenberg number	W e	0.1, 0.5, 0.9	F s = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , P r = 0.9 , ϵ 1 = 0.1	Figure 2
2	Forchheimer number	F s	0.1, 0.5, 0.9	W e = 0.1 , D a = 7 , N = 0.1 , M = 0.1 , P r = 0.9 , ϵ 1 = 0.1	Figure 3
3	Buoyancy ratio	N	0.1, 0.5, 0.9	W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , P r = 0.9 , ϵ 1 = 0.1	Figure 4
4	Thermal conductivity parameter	ϵ 1	0.1, 0.5, 0.9	W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , P r = 0.9 , N = 0.1	Figure 5
5	Reaction rate parameter	k c	0.1, 0.3, 0.9	W e = 0.1 , D a = 5 , F s = 0.1 , M = 0.1 , P r = 0.9 , S c = 0.9 , ϵ 1 = 0.1	Figure 6
6	Prandtl and ϵ 1	P r , ϵ 1	0.7–1.1	Other parameters are fixed as above	Figure 7
7	Schmidt and the reaction rate	S c , k c	0.7–1.1	Other parameters are fixed as above	Figure 8

5.9 Computational cost comparison

As shown in Table 3, for a fixed grid size of 50 × 500 , the proposed scheme consumes 30–40% less runtime than the classical second-order Runge–Kutta method while achieving superior accuracy (lower L 2 -norm error). Moreover, due to the higher order of spatial accuracy, fewer grid points are required to attain the same resolution, reducing computational cost in long-time simulations.

Table 3

Comparison of the proposed scheme and second-order Runge–Kutta method for Stokes’ first problem ( N x ( grid points ) = 50, final time t f = 1 ; L 2 error, and runtime for various time steps)

Time steps N t	Scheme type	Spatial scheme	L 2 -error	Runtime (s)
250	Proposed scheme	Sixth-order compact	0.0747	0.78
	Second-order Runge–Kutta	Sixth-order compact	0.0770	1.13
300	Proposed scheme	Sixth-order compact	0.0732	0.85
	Second-order Runge–Kutta	Sixth-order compact	0.0751	1.26
350	Proposed scheme	Sixth-order compact	0.0730	0.93
	Second-order Runge–Kutta	Sixth-order compact	0.0746	1.34
400	Proposed scheme	Sixth-order compact	0.0734	1.06
	Second-order Runge–Kutta	Sixth-order compact	0.0748	1.52
450	Proposed scheme	Sixth-order compact	0.0743	1.22
	Second-order Runge–Kutta	Sixth-order compact	0.0755	1.68
500	Proposed scheme	Sixth-order compact	0.0753	1.39
	Second-order Runge–Kutta	Sixth-order compact	0.0763	1.88

6 Conclusion

This research developed and assessed a robust numerical system to model the Darcy–Forchheimer flow of a Williamson fluid whose thermal conductivity varies with temperature. The proposed scheme combines a second-order time-accurate computational technique with a sixth-order compact spatial discretization. The system is conditionally stable, as shown by von Neumann stability analysis or the Fourier series. This provides a solid foundation for its efficiency in solving complex time-dependent partial differential equations. An explicit scheme combined the existing Runge–Kutta scheme and a modified exponential integrator. The scheme comprises two stages: the predictor and corrector stages. The predictor phase employs a modified exponential integrator, while the corrector method utilizes the second order of the established Runge–Kutta scheme. The governing equations for heat and mass transfer of Williamson fluids have been transformed to dimensionless partial differential equations. A modification of the proposed scheme was suggested compared with one of the existing schemes. The numerical results show that the suggested scheme outperforms the second-order Runge–Kutta approach regarding accuracy and error rate. The combination of sophisticated time and spatial discretization techniques, which successfully deal with the mathematical model’s non-linearities and complex relationships, is responsible for this advantage. The proposed scheme was conditionally stable, and its convergence was mathematically proven. Quantitative validation against the classical second-order Runge–Kutta method revealed that the proposed method reduces error norms by up to 18% and maintains accuracy even for small time-step sizes. Incorporating a change into the suggested system also shows how flexible it is and how much room there is for improvement in different computational contexts. The scheme successfully captured physical trends, as follows:

The scheme was better than the existing second-order Runge–Kutta scheme for finding the error norm for Stokes’ first problem.
The velocity profile diminished as the Weissenberg parameter increased.
The increase in Forchheimer numbers produced decay in the velocity profile.
The temperature profile had dual behaviour by uplifting small dimensionless parameters in variable thermal conductivity.

This research addresses a significant gap in the numerical modelling of non-Newtonian fluid flows that exhibit variations in thermal properties. This work provides valuable insights for researchers and engineers engaged in fluid flow in porous media, heat transfer processes, and associated industrial applications.

In the future, researchers may investigate how to make the suggested method work with three-dimensional flow models that include more complicated things like changing porosity, chemical reactions, or multiphase flow. To improve stability and performance for larger time-step applications, the creation of fully implicit or semi-implicit schemes can also be investigated. These improvements would make the suggested computational framework even more valuable and essential for solving fluid flow problems in the real world. Various engineering and industrial applications are relevant to the proposed model and numerical framework. The Darcy–Forchheimer–Williamson formulation with variable thermal conductivity adequately describes physical phenomena in increased oil recovery, polymer extrusion processes, cooling of electronic devices through porous heat sinks, and blood flow through biological tissues under magnetic fields. The model is relevant to solar collectors, electromagnetic pumps, and geothermal energy systems, where thermal gradients and electromagnetic fields affect flow behaviour since they include mixed convection and MHD effects. The high-accuracy numerical scheme created in this work supports improved design and optimization in real applications by allowing quick simulation of such complicated systems.

Acknowledgments

The author would like to acknowledge the support of Prince Sultan University for paying the article processing charges of this publication.

Funding information: This research was supported by Prince Sultan University for paying the article processing charges of this publication.
Author contributions: The author has accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The author states no conflict of interest.
Data availability statement: Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

References

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

[2] Nadeem S, Akram S. Influence of inclined magnetic field on peristaltic flow of a Williamson fluid model in an inclined symmetric or asymmetric channel. Math Comput Model. 2010;52:107–19.10.1016/j.mcm.2010.02.001Search in Google Scholar

[3] Nadeem S, Hussain ST, Lee C. Flow of a Williamson fluid over a stretching sheet. Braz J Chem Eng. 2013;30(3):619–25.10.1590/S0104-66322013000300019Search in Google Scholar

[4] Akbar NS, Nadeem S, Lee C, Khan ZH, Haq RU. Numerical study of Williamson nanofluid flow in an asymmetric channel. Results Phys. 2013;3:161–6.10.1016/j.rinp.2013.08.005Search in Google Scholar

[5] Eldabe NT, Elogail MA, Elshaboury SM, Hasan AA. Hall effects on the peristaltic transport of Williamson fluid through a porous medium with heat and mass transfer. Appl Math Model. 2016;40:315–28.10.1016/j.apm.2015.04.043Search in Google Scholar

[6] Krishnamurthy MR, Prasannakumara BC, Gireesha BJ, Gorla RSR. Effect of chemical reaction on MHD boundary layer flow and melting heat transfer of Williamson nanofluid in porous medium. Eng Sci Technol Int J. 2016;19:53–61.10.1016/j.jestch.2015.06.010Search in Google Scholar

[7] Bhatti MM, Rashidi MM. Effects of thermo-diffusion and thermal radiation on Williamson nanofluid over a porous shrinking/stretching sheet. J Mol Liq. 2016;221:567–73.10.1016/j.molliq.2016.05.049Search in Google Scholar

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

[9] Kumaran G, Sandeep N. Thermophoresis and Brownian moment effects on parabolic flow of MHD Casson and Williamson fluids with cross diffusion. J Mol Liq. 2017;233:262–9.10.1016/j.molliq.2017.03.031Search in Google Scholar

[10] Hayat T, Saleem A, Tanveer A, Alsaadi F. Numerical study for MHD peristaltic flow of Williamson nanofluid in an endoscope with partial slip and wall properties. Int J Heat Mass Transf. 2017;114:1181–7.10.1016/j.ijheatmasstransfer.2017.06.066Search in Google Scholar

[11] Hamid A, Hashim I, Khan M. Numerical simulation for heat transfer performance in unsteady flow of Williamson fluid driven by a wedge-geometry. Results Phys. 2018;9:479–85.10.1016/j.rinp.2018.01.025Search in Google Scholar

[12] Hashim I, Khan M, Hamid A. Numerical investigation on time-dependent flow of Williamson nanofluid along with heat and mass transfer characteristics past a wedge geometry. Int J Heat Mass Transf. 2018;118:480–91.Search in Google Scholar

[13] Hashim I, Hamid A, Khan M. Unsteady mixed convective flow of Williamson nanofluid with heat transfer in the presence of variable thermal conductivity and magnetic field. J Mol Liq. 2018;260:436–46.10.1016/j.molliq.2018.03.079Search in Google Scholar

[14] Bahiraei M, Mazaheri N, Alighardashi M. Development of chaotic advection in laminar flow of a non-Newtonian nanofluid: A novel application for efficient use of energy. Appl Therm Eng. 2017;124:1213–23.10.1016/j.applthermaleng.2017.06.106Search in Google Scholar

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

[16] Mabood F, Khan WA, Makinde OD. Hydromagnetic flow of a variable viscosity nanofluid in a rotating permeable channel with Hall effects. J Eng Thermophys. 2017;2:553–66.10.1134/S1810232817040105Search in Google Scholar

[17] Zafar SS, Zaib A, Faizan M, Shah NA, Ali F, Yook SJ. Bioconvection flow of Prandtl nanomaterial due to stretched cylinder enclosed through Darcy–Forchheimer flow with triple stratification. Alex Eng J. 2025;116:188–201.10.1016/j.aej.2024.12.062Search in Google Scholar

[18] Reddy MV, Nath JM, Ali F, Das TK, Khan U, Garayev M. Implementation of homotopy analysis method for entropy-optimized two-phase nanofluid flow in a bioconvective non-Newtonian model with thermal radiation. J Radiat Res Appl Sci. 2025;18(1):101.10.1016/j.jrras.2024.101218Search in Google Scholar

[19] Rehman KU, Shatanawi W, Abodayeh K. Thermophysical aspects of magnetized Williamson fluid flow subject to both porous and non-porous surfaces: A Lie symmetry analysis. Case Stud Therm Eng. 2021;28:101688.10.1016/j.csite.2021.101688Search in Google Scholar

[20] Merkin JH. A model for isothermal homogeneous–heterogeneous reactions in boundary layer flow. Math Comput Model. 1996;24:125–36.10.1016/0895-7177(96)00145-8Search in Google Scholar

[21] Chaudhary MA, Merkin JH. A simple isothermal model for homogeneous–heterogeneous reactions in boundary layer flow: I. Equal diffusivities. Fluid Dyn Res. 1995;16:311–3.10.1016/0169-5983(95)00015-6Search in Google Scholar

[22] Ziabakhsh Z, Ganji DD, Bararnia H, Babazadeh H. Analytical solution of flow and diffusion of chemically reactive species over a nonlinearly stretching sheet immersed in a porous medium. J Taiwan Inst Chem Eng. 2010;41:22–8.10.1016/j.jtice.2009.04.011Search in Google Scholar

[23] Khan WA, Pop I. Flow near the two-dimensional stagnation-point on an infinite permeable wall with a homogeneous–heterogeneous reaction. Commun Nonlinear Sci Numer Simul. 2010;15:3435–43.10.1016/j.cnsns.2009.12.022Search in Google Scholar

[24] Kameswaran PK, Shaw S, Sibanda P, Murthy PVSN. Homogeneous–heterogeneous reactions in a nanofluid flow due to a porous stretching sheet. Int J Heat Mass Transf. 2013;57:465–72.10.1016/j.ijheatmasstransfer.2012.10.047Search in Google Scholar

[25] Hayat T, Hussain Z, Alsaedi A, Asghar S. Carbon nanotubes effects in the stagnation point flow towards a non-linear stretching sheet with variable thickness. Adv Powder Technol. 2016;27:1677–88.10.1016/j.apt.2016.06.001Search in Google Scholar

[26] Chen J, Liu B, Gao X, Yan L, Xu D. Effects of heterogeneous–homogeneous interaction on the homogeneous ignition in hydrogen-fueled catalytic microreactors. Int J Hydrog Energy. 2016;41:11441–54.10.1016/j.ijhydene.2016.05.022Search in Google Scholar

[27] Hashim I, Khan M. On Cattaneo–Christov heat flux model for Carreau fluid flow over a slendering sheet. Results Phys. 2017;79:310–9.10.1016/j.rinp.2016.12.031Search in Google Scholar

[28] Khan MI, Waqas M, Hayat T, Khan MI, Alsaedi A. Numerical simulation of non-linear thermal radiation and homogeneous–heterogeneous reactions in convective flow by a variable thicked surface. J Mol Liq. 2017;46:259–67.10.1016/j.molliq.2017.09.075Search in Google Scholar

[29] Khan M, Ahmad L, Khan WA, Alshomrani AS, Alzahrani AK, Alghamdi MS. A 3D Sisko fluid flow with Cattaneo–Christov heat flux model and heterogeneous–homogeneous reactions: A numerical study. J Mol Liq. 2017;238:19–26.10.1016/j.molliq.2017.04.059Search in Google Scholar

[30] Khan RA, Ullah H, Raja MAZ, Khan MAR, Islam S, Shoaib M. Heat transfer between two porous parallel plates of steady nanofluids with Brownian and thermophoretic effects: A new stochastic numerical approach. Int Commun Heat Mass Transf. 2021;126:105436.10.1016/j.icheatmasstransfer.2021.105436Search in Google Scholar

[31] Ullah H, Khan I, AlSalman H, Islam S, Raja MAZ, Shoaib M, et al. Levenberg–Marquardt backpropagation for numerical treatment of micropolar flow in a porous channel with mass injection. Complexity. 2021;2021:5337589.10.1155/2021/5337589Search in Google Scholar

[32] Shoaib M, Khan RA, Ullah H, Nisar KS, Raja MAZ, Islam S, et al. Heat transfer impacts on Maxwell nanofluid flow over a vertical moving surface with MHD using stochastic numerical technique via artificial neural networks. Coatings. 2021;11(12):1483.10.3390/coatings11121483Search in Google Scholar

[33] Rehman KU, Shatanawi W, Alharbi WG, Shatnawi TA. AI-neural networking analysis (NNA) of thermally slip magnetized Williamson (TSMW) fluid flow with heat source. Case Stud Therm Eng. 2024;56:104248.10.1016/j.csite.2024.104248Search in Google Scholar

[34] 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-9Search in Google Scholar PubMed PubMed Central

[35] Williamson RV. The flow of pseudoplastic materials. Ind Eng Chem Res. 1929;21(11):1108–11.10.1021/ie50239a035Search in Google Scholar

[36] Hashim M, Khan M, Hamid A. Numerical investigation on time-dependent flow of Williamson nanofluid along with heat and mass transfer characteristics past a wedge geometry. Int J Heat Mass Transf. 2018;118:480–91.10.1016/j.ijheatmasstransfer.2017.10.126Search in Google Scholar

[37] Khan M, Salahuddin T, Malik MY, Mallawi FO. Change in viscosity of Williamson nanofluid flow due to thermal and solutal stratification. Int J Heat Mass Transf. 2018;126:941–8.10.1016/j.ijheatmasstransfer.2018.05.074Search in Google Scholar

[38] Salahuddin T. Modelling unsteady waveform Williamson fluid flow near a permeable radioactive surface. Int Commun Heat Mass Transf. 2020;117:104764.10.1016/j.icheatmasstransfer.2020.104764Search in Google Scholar

[39] Bilal M, Sagheer M, Hussain S, Mehmood Y. MHD stagnation point flow of Williamson fluid over a stretching cylinder with variable thermal conductivity and homogeneous/heterogeneous reaction. Commun Theor Phys. 2017;67(6):688–96.10.1088/0253-6102/67/6/688Search in Google Scholar

[40] Subbarayudu K, Suneetha S, Reddy PBA. The assessment of time dependent flow of Williamson fluid with radiative blood flow against a wedge. Propul Power Res. 2020;9(1):87–99.10.1016/j.jppr.2019.07.001Search in Google Scholar

[41] Ganesh Kumar K, Rudraswamy NG, Gireesha BJ, Krishnamurthy MR. Influence of non-linear thermal radiation and viscous dissipation on three-dimensional flow of Jeffrey nanofluid over a stretching sheet in the presence of Joule heating. Nonlinear Eng. 2017;6(3):207–19.10.1515/nleng-2017-0014Search in Google Scholar

[42] Usman MI, Khan F, Shah SU, Khan A, Ghaffari Y, Chu YM. Heat and mass transfer analysis for bioconvective flow of Eyring Powell nanofluid over a Riga surface with non-linear thermal features. Numer Methods Partial Differ Equ. 2020;36(9):1803–21.Search in Google Scholar

[43] Shehzad SA, Hayat T, Alsaedi A, Obid MA. Non-linear thermal radiation in three-dimensional flow of Jeffrey nanofluid: A model for solar energy. Appl Math Comput. 2014;248:273–86.10.1016/j.amc.2014.09.091Search in Google Scholar

[44] Hayat T, Imtiaz M, Alsaedi A, Kutbi MA. MHD three-dimensional flow of nanofluid with velocity slip and non-linear thermal radiation. J Magn Magn Mater. 2015;396:31–7.10.1016/j.jmmm.2015.07.091Search in Google Scholar

[45] Usman MI, Khan SU, Khan A, Ghaffari Y, Chu YM, Farooq S. A higher order slip flow of generalized micropolar nanofluid with applications of motile microorganisms, non-linear thermal radiation and activation energy. Int J Mod Phys B. 2021;35(7):2150095.10.1142/S0217979221500958Search in Google Scholar

[46] Animasaun IL, Raju CSK, Sandeep N. Unequal diffusivities case of homogeneous-heterogeneous reactions within viscoelastic fluid flow in the presence of induced magnetic field and non-linear thermal radiation. Alex Eng J. 2016;55(2):1595–606.10.1016/j.aej.2016.01.018Search in Google Scholar

[47] Azam M. Effects of Cattaneo-Christov heat flux and non-linear thermal radiation on MHD Maxwell nanofluid with Arrhenius activation energy. Case Stud Therm Eng. 2022;34:102048.10.1016/j.csite.2022.102048Search in Google Scholar

Received: 2025-01-24

Revised: 2025-05-15

Accepted: 2025-07-01

Published Online: 2025-07-25

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

Articles in the same Issue

https://doi.org/10.1515/phys-2025-0190

Keywords for this article

computational scheme; stability analysis; convergence behaviour; Darcy–Forchheimer flow; Williamson fluid; heat and mass transfer; porous media; mixed convection; numerical comparison

Creative Commons

BY 4.0