In recent years, the investigation of non-Newtonian fluids has gained considerable interest owing to its complex behaviour, which is crucial in numerous engineering and commercial applications, including petroleum extraction, polymer processing, and biological systems. The Williamson fluid is a non-Newtonian fluid distinguished by its shear-thinning characteristics. The flow characteristics of Williamson fluids, particularly when affected by stochastic factors, pose distinct challenges in modelling and computation because of the nonlinear structure of their governing equations. Numerical strategies are essential for achieving approximate solutions when analytical methods are impractical.

We formulate a resilient two-stage predictor-corrector methodology to tackle the flow’s stochastic characteristics. The predictor phase calculates the solution at an intermediate temporal level, yielding a preliminary estimate, while the corrector phase enhances this solution at the following temporal level. This method facilitates precise time-stepping for the stochastic differential equations that dictate the flow of Williamson fluid, rendering it appropriate for managing the system’s intrinsic uncertainties and unpredictable variations.

Here are the main points of this work:

1. We developed a two-stage numerical scheme for solving deterministic and stochastic partial differential equations with second-order temporal and sixth-order spatial accuracy.

Non-Newtonian fluids exhibit distinct rheological characteristics, defined by their shear-thinning or shear-thickening behaviour. The characteristics mentioned above confer benefits in domains such as fluid transport, mixing, and heat transfer, thereby enhancing the efficiency and performance of various processes. Consequently, these fluids play an essential role in shaping and advancing the domain of modern engineering and industrial methodologies. Recently, a growing body of research has focused on a range of non-Newtonian models, encompassing the viscoelastic system [1], Maxwell system [2], Casson system [3], power-law system [4], Carreau system [5], and Eyring-Powell system [6].

The Williamson model, a significant element in emulsion sheet manufacture, exemplifies one category of such models. This notion is applied in various contexts, including plasma flow, photographic film production, and the understanding of hemodynamics. The shear-thinning properties of the Williamson [7] model in non-Newtonian fluids are well-known. He stressed how this model is crucial for differentiating plastic flow from viscous flow. In addition, he found that this model is just as important in bioengineering, especially for evaluating hemodialysis applications and mass and heat transport in blood arteries. Due to its significance, numerous scholars have dedicated their time and resources to investigating the dynamics of Williamson fluid flow across various conditions, demonstrating the research community’s deep interest and commitment to comprehensively understanding this fluid type. Owing to their importance, these fluids are essential in numerous vital technical and industrial domains. The homotopy analysis method can be used for quantitative analysis [8], in stagnation point flow [9], chemical reactions on a stretching cylinder [10], in porous media [11], with the effects of viscous dissipation and slip velocity [12], with thermal radiation and chemical reaction influences [13], under multiple slip boundary conditions [14], and in nanofluid scenarios with magnetic field effects [15].

According to Williamson, the Williamson fluid model is one of the simplest non-Newtonian models that attempts to capture the behaviour of viscoelastic shear-thinning [16]. A chemical process was used by Krishnamurthy et al. [17] to demonstrate the flow of thermally radiative Williamson fluid on a stretched sheet. Their findings demonstrated a decrease in fluid temperature attributable to the existence of the Williamson parameter. Khan et al. [18] illustrated the effects of slip flow of Williamson nanofluid within a porous media. The surface drag force diminishes as the Williamson fluid parameter increases. Hayat et al. [19] examined a Williamson fluid’s 2D unsteady radiative flow on a porous stretched surface. As the Williamson parameter increases, the fluid velocity decreases. In their investigation of Williamson fluid flow across a stretching sheet, Nadeem et al. [20] discovered that the skin friction coefficient decreases as the Williamson parameter increases. Apply the Keller box method to investigate the MHD flow of Williamson fluid across a stretching sheet, as demonstrated by Salahuddin et al. [21]. Their results suggest that the Williamson fluid parameter decreases fluid velocity. Refs. [22, 23] contain only a few substantial analyses of this subject matter.

Many industrial processes include fluids passing through porous media. Its few uses include drying wood, storing nuclear waste, processing food, refining oil, facilitating drainage, and watering crops. Darcy’s principle analyses flow behaviour under low porosity and small velocity conditions. The Darcy principle didn't work when the Reynolds number amount was more than 1. Forchheimer [24] included the square velocity element in the momentum equation to get around this restriction. Then, for operations involving higher Reynolds numbers, this is referred to as the Forchheimer number. The Darcy–Forchheimer flow of a viscous fluid across a plate was examined by Mukhopadhyay et al. [25] using numerical analysis. The researchers found that lowering the permeability parameter made the fluid colder. The 3D Williamson nanomaterial flow over a Darcy–Forchheimer porous media is shown by Hayat et al. [26]. As the Forchheimer number increases, the surface shear stress decreases. Khan et al. [27] visualized the Darcy–Forchheimer flow of a viscous fluid undergoing heterogeneous-homogeneous chemical processes. The availability of the Darcy number causes a noticeable reduction in the fluid speed, as their data demonstrate. An analysis of the hybrid nanofluid’s Darcy–Forchheimer and slip flow on a revolving disc is carried out by Haider et al. [28]. They demonstrated that a greater estimate of Forchheimer enhances the fluid temperature. Sadiq et al. [29] identified a steady 3D Darcy–Forchheimer flow of carbon nanotubes on a spinning disc in Refs. [30–33], you can find a compilation of significant research on these topics.

A viscous fluid that cannot be compressed and is incompressible is described by the Navier-Stokes equation. Despite the notion that the Navier-Stokes equation is intrinsically deterministic, the fact that certain solutions behave randomly may lead to fresh insights regarding the beginning of turbulence. However, from a purely mathematical point of view, whether solutions to the equation presented above exist and whether they are smooth is still largely unsolved. In the following paragraphs, we will discuss several approaches to analyzing the solutions to the above equation, including a random component.

One famous example of a deterministic equation that conceals unpredictability is the famous Schrödinger equation. Any quantum physics experiment will show some degree of blatant randomness, even though it is not based on mathematical probability theory. Notwithstanding this inexplicable dilemma, Quantum Theory has yielded an impressive suite of control tools, notably the inversion of solutions to classical equations of motion that are (in principle) differentiable and instead display uncertainty. Simulating “dry water” (Von Neumann) only adds another layer of complexity to the already complicated Euler equation.

There seems to be no proof currently that it does not produce singularities [34]. There are other ways to generate randomness, although errors in the initial conditions could also be a source of uncertainty. This calls for applying statistical approaches, such as tracking the evolution of a probability measure in conjunction with the pertinent physical beginning data. See, for instance [35, 36], for examples of how this fits into the statistical approach to turbulence that began in the 19th century. A variety of stochastic diffusion-based Langevin dynamics may account for both equilibrium and non-equilibrium dynamics, as well as Kraichnan’s model in turbulent advection [37]. Nevertheless, the chosen numerical model may include uncertainty (see [38, 39] for further information on this subject in climate modelling).

A popular method for dealing with stochasticity is introducing stochastic partial differential equations by applying random forces to the Navier-Stokes equation. After the revolutionary mathematical work, a deluge of material covers the subject [40]. Even though turbulence usually introduces stochasticity at the Eulerian level, models of the Langevin type that include smooth Lagrangian trajectories and stochastic velocities have been evaluated [41]. D.D. Holm more recently used Lie transfer to introduce stochastic advection [42]. Stochastic partial differential equations describe the resultant motion, and the approach is also Eulerian.

This space cannot accommodate any stochastic methods related to fluid dynamics. A limited number of subjects and their corresponding references have been chosen. It is important to note that several fascinating formulas exist for probabilistic representations. This is a recognized practice in stochastic analysis. A potential application lies in fluid dynamics, where it may be utilized to depict the anticipated values of functionals of stochastic processes as solutions to partial differential equations. Three distinct methodologies are presented: one employs a probabilistic model of the vorticity field, another utilizes branching processes alongside the Fourier transform, and the third implements Lagrangian diffusion processes. For additional information, go to [43–45] accordingly.

This study offers two numerical schemes that can be used to discretize partial differential equations. The scheme is constructed using Taylor series expansion and then implemented to solve the flow problem. The non-Newtonian, incompressible, laminar fluid flow is considered over the flat plate. The governing equations are reduced to dimensionless partial ones rather than ordinary ones. Then, these dimensionless equations are solved using the proposed scheme.

The proposed scheme is applicable for solving stochastic differential equations; however, a scheme for the deterministic model must be constructed before its construction. The scheme is a bifurcated predictor-corrector framework. The predictor scheme determines the answer at a specified time level, whereas the corrector scheme identifies the solution at the subsequent time level. To propose a scheme, consider the following differential equations. ∂ v ∂ t = γ ∂ 2 v ∂ y 2 + f v (1) where f v represents the nonlinear term in the system, this framework ensures that nonlinear and stochastic effects are managed systematically, avoiding excessive numerical diffusion and providing better stability.

The Von Neumann stability analysis, or Fourier series analysis, is used to assess the stability of finite difference schemes. Using this stability analysis, the difference equation is transformed into trigonometric equations. Further stability conditions are imposed into the trigonometric equations. The analysis provides the exact condition for linear partial differential equations and estimates the actual stability condition for nonlinear partial differential equations. To apply this analysis, consider the following transformations C e i I ψ = β 1 e i + 1 I ψ + e i I ψ + β 1 e i − 1 I ψ (16) D e i I ψ = b ∘ e i + 1 I ψ − 2 e i I ψ + e i − 1 I ψ Δ y 2 + b 1 e i + 2 I ψ − 2 e i I ψ + e i − 2 I ψ 4 Δ y 2 (17) Where I = − 1 is an imaginary number.

Applying transformations (Equation 16) and (Equation 17) into the predictor stage of the proposed scheme (Equation 11) yields v ¯ i n + 1 = v i n e 0.05 Δ t + e 0.05 Δ t − 1 0.05 γ   4 b ∘ cos ⁡ ψ − 1 + 2 b 1 cos ⁡ ψ − 1 2 Δ y 2 2 β 1 ⁡ cos ⁡ ψ + 1 − 0.05 v i n (18)

Re-write Equation 18 as v ¯ i n + 1 = δ v i n (19) Where δ = e 0.05 Δ t + e 0.05 Δ t − 1 0.05 γ   4 b ∘ cos ⁡ ψ − 1 + 2 b 1 cos ⁡ ψ − 1 2 Δ y 2 2 β 1 ⁡ cos ⁡ ψ + 1 − 0.05

Incorporating transformation (Equations 16, 17) into the second corrector stage of the scheme (Equation 15) gives using f = 0 . v i n + 1 = a v i n + b v ¯ i n + 1 + c e Δ t − 1 γ 4 b ∘ cos ⁡ ψ − 1 + 2 b 1 cos ⁡ ψ − 1 2 Δ y 2 2 β 1 ⁡ cos ⁡ ψ + 1 v ¯ i n + 1 + σ v i n Δ W (20)

By using the first stage (Equation 19) in the second stage (Equation 20) it yields v i n + 1 = a v i n + b δ v i n + c e Δ t − 1 γ 4 b ∘ cos ⁡ ψ − 1 + 2 b 1 cos ⁡ ψ − 1 2 Δ y 2 2 β 1 ⁡ cos ⁡ ψ + 1 δ v i n + σ v i n Δ W (21)

Re-write Equation 21 as v i n + 1 = δ 1 v i n + σ v i n Δ W (22) Where δ 1 = a + b δ + c e Δ t − 1 γ 4 b ∘ cos ⁡ ψ − 1 + 2 b 1 cos ⁡ ψ − 1 2 Δ y 2 2 β 1 ⁡ cos ⁡ ψ + 1 δ in Equation 22

The amplification factor can be expressed as v i n + 1 v i n = δ 1 + σ Δ W (23)

By using the expected value of the square of the amplification factor, Equation 23 is written as E v i n + 1 v i n 2 ≤ 2 E δ 1 2 + 2 σ 2 E Δ W 2 (24)

If δ 1 2 < 1 and let, 2 σ 2 = λ then inequality (Equation 24) is stated as E v i n + 1 v i n 2 ≤ 1 + λ Δ t (25)

The scheme will remain stable in the mean square sense if it meets conditions (Equation 25); otherwise, the solution will be unstable.

Theorem 1

The proposed time and compact scheme in space are consistent in the mean square sense.

Proof: Let P be the smooth function L P i n = P n + 1 Δ t , i Δ x − P n Δ t , i Δ x − γ ∫ n Δ t n + 1 Δ t P y y s , i Δ x d S − σ ∫ n Δ t n + 1 Δ t P s , i Δ x d W s (26) L i n P = P n + 1 Δ t , i Δ x − P n Δ t , i Δ x − b e 0.05 Δ t − 1 0.05 C − 1 D P n Δ t , i Δ x + c e Δ t − 1 C − 1 D P ¯ n + 1 Δ t , i Δ x − σ P n Δ t , i Δ x W n + 1 Δ t − W n Δ t (27) Where P ¯ n + 1 Δ t , i Δ x = P n Δ t , i Δ x e 0.05 Δ t + e 0.05 Δ t − 1 0.05 { C − 1 D P n Δ t , i Δ x − 0.05 P ( n Δ t , i Δ x )

By subtracting Equation 27 from Equation 26 and then applying the absolute value of the square of difference it yields E L P i n − L i n P 2 = E − γ ∫ n Δ t n + 1 Δ t P y y s , i Δ x d S − σ ∫ n Δ t n + 1 Δ t P s , i Δ x d W s + b e 0.05 Δ t − 1 0.05 C − 1 D P n Δ t , i Δ x + c e Δ t − 1 C − 1 D P ¯ n + 1 Δ t , i Δ x + σ P n Δ t , i Δ x W n + 1 Δ t 2 (28)

Re-write Equation 28 as E L P i n − L i n P 2 ≤ 2 γ E ∫ n Δ t n + 1 Δ t P y y s , i Δ x d S − b e 0.05 Δ t − 1 0.05 C − 1 D P n Δ t , i Δ x + c e Δ t − 1 C − 1 D P ¯ n + 1 Δ t , i Δ x 2 + 2 σ 2 E ∫ n Δ t n + 1 Δ t P s , i Δ x d W s − P n Δ t , i Δ x W n + 1 Δ t − W n Δ t 2 (29)

Now consider the inequality E ∫ t ∘ t f s , W d W s 2 m ≤ t − t ∘ m − 1 m 2 m − 1 m ∫ t ∘ t E f s , W 2 m d s (30)

The use of inequality (Equation 30) in inequality (Equation 29) yields E L P i n − L i n P 2 ≤ 2 γ E ∫ n Δ t n + 1 Δ t P y y s , i Δ x d S − b e 0.05 Δ t − 1 0.05 C − 1 D P n Δ t , i Δ x + c e Δ t − 1 C − 1 D P ¯ n + 1 Δ t , i Δ x 2 + 2 σ 2 Δ t ∫ n Δ t n + 1 Δ t E P s , i Δ x − P n Δ t , i Δ x 2 dS (31)

Upon applying the l i m i t in Equation 31 as Δ x → 0 , Δ t → 0 , t = n + 1 Δ t , n Δ t , i Δ x → t , x , then the mean square error tends to zero.i.e. E L P i n − L i n P 2 → 0 (32)

Thus, the proposed scheme in time and compact spatial discretization is consistent in the mean square sense proved by Equation 32.

Here, we provide a comparative summary highlighting the stability criteria of our proposed scheme versus alternative methods. Below, we outline a summary of Table 1 comparing the stability characteristics of the proposed method against existing numerical schemes.

The Von Neumann stability analysis is applied to derive the stability condition. The stochastic nature of the problem requires the scheme to be stable in the mean-square sense. The derived stability condition is: E v i n + 1 v i n 2 ≤ 1 + λ Δ t . This ensures the numerical scheme remains bounded over time, preventing error amplification.

Table 1 effectively compares the stability criteria, stochastic effects, accuracy, and efficiency of the proposed method versus alternative approaches. The proposed scheme ensures stability in stochastic problems and offers superior accuracy and computational efficiency.

Examine a laminar, unstable, incompressible Williamson fluid above a fixed plate. The fluid flow is propelled by temperature and concentration gradients. Let x * -axis is taken perpendicular along the plate and y * - axis is taken horizontally. Suppose the ambient temperature and concentration are less than those at the plate. At t * = 0 the temperature and concentration are represented by T ∞ and C ∞ respectively. After t * > 0 both temperature and concentration become T ∞ and C ∞ plus some periodic functions. Figure 1 illustrates the physical and mathematical configuration of the problem involving non-Newtonian Williamson fluid flow over a stationary vertical sheet under the influence of Darcy–Forchheimer porous medium effects, thermal diffusion, and solutal diffusion. The x * -axis represents the horizontal direction along the sheet. The y * -axis represents the normal direction (perpendicular) to the sheet. The coordinate system is essential in defining velocity components, boundary layer thicknesses, and external forces acting on the flow. The figure depicts three different boundary layers: The momentum boundary layer (Black Curve) Represents the region where the velocity of the fluid gradually transitions from zero (no slip at the wall) to the free-stream velocity. Thermal boundary layer (Orange Curve): Shows the temperature distribution within the fluid, where heat transfer occurs from the wall to the fluid. Concentration boundary layer (Purple Curve): Indicates the distribution of solute concentration in the fluid, where mass transfer occurs due to concentration gradients. Gravity ( g : Acts downward, influencing the buoyancy-driven (natural convection) effects in the fluid flow. Magnetic Field B ∘ : Represents an applied transverse magnetic field, which affects the flow characteristics due to the magnetohydrodynamics (MHD) effect. Porous Medium Effects: The presence of Darcy–Forchheimer drag forces affects the momentum balance, introducing resistance to the fluid flow. At the wall (y* = 0): The velocity u * is zero due to the no-slip condition, and the temperature T and concentration C are prescribed. The expressions: u * = 0 , T = T w + ϵ 1 T w − T ∞ cos w t , C = C w + ϵ 1 C w − C ∞ cos w t indicate that temperature and concentration vary periodically over time due to oscillations in wall properties. Far from the wall (as y * → ∞ ): The velocity approaches zero, and temperature and concentration reach their respective free-stream values T ∞ and C ∞ . Considering the effect of the Darcy Forchheimer model, the governing equations of the flow can be expressed as [46] ∂ u * ∂ t * = ν ∂ 2 u * ∂ y * 2 + 2 ν Γ ∂ 2 u * ∂ y * 2 ∂ u * ∂ y * − b ¯ u * 2 + g   β T T − T ∞ + g β c C − C ∞ + σ B ∘ 2 ρ u * − ν k p u * (33) ∂ T ∂ t * = α ∂ 2 T ∂ y * 2 + q ‴ (34) ∂ C ∂ t *   = D ∂ 2 C ∂ y * 2 − k c C − C ∞ (35)

Here is the explanation of each equation used in problem formulation:

Momentum (Velocity) Equation 33: ∂ u * ∂ t * : represents the time-dependent change in velocity. ν ∂ 2 u * ∂ y * 2 : represents the viscous diffusion term (Laplacian of velocity). 2 ν Γ ∂ 2 u * ∂ y * 2 ∂ u * ∂ y * : Accounts for the non-Newtonian Williamson fluid effects, where Γ is the Williamson parameter. b ¯ u * 2 : represents the Forchheimer drag force, a correction to the Darcy resistance term. g   β T T − T ∞ : represents buoyancy forces due to thermal variations. g β c C − C ∞ : represents buoyancy forces due to concentration variations. σ B ∘ 2 ρ u * : represents the magnetohydrodynamic (MHD) Lorentz force, which opposes the flow. − ν k p u * : represents Darcy’s resistance in the porous medium.

Energy (Heat Transfer) Equation 34: ∂ T ∂ t * : represents the time-dependent change in temperature. α ∂ 2 T ∂ y * 2 : represents heat conduction (thermal diffusion). q ‴ = k u R ρ x * ν c p A * u T w − T ∞ + B * T − T ∞ : represents the internal heat generation where A * and B * are dimensionless coefficients of space and temperature-dependent terms, respectively.

Concentration (Mass Diffusion) Equation 35: ∂ C ∂ t *   : represents the time-dependent change in concentration. D ∂ 2 C ∂ y * 2 : represents mass diffusion. − k c C − C ∞ : Represents a first-order chemical reaction, where reactants are consumed over time.

Subject to initial and boundary conditions u * = 0 , T = T ∞ , C = C ∞   w h e n   t * = 0 u * = 0 , T = T ∞ + ϵ 1 T w − T ∞ cos ⁡ w * t * , C = C ∞ + ϵ 1 C w − C ∞ cos ⁡ w * t *   w h e n   y * = 0 u * → 0 , T → T ∞ , C → C ∞   w h e n   y * → ∞ (36)

Let B ∘ is the strength of the magnetic field applied transversely to the plate. Now, using the transformations (Equation 36). By introducing non-dimensional variables, the governing equations are transformed into a more convenient form for numerical analysis:

We present a detailed simulation study with the following objectives:

This research introduces an innovative computational method for addressing deterministic and stochastic partial differential equations, specifically emphasizing the stochastic Darcy–Forchheimer flow of Williamson fluid over a stationary surface. The suggested two-stage approach integrates a modified time integrator with a second-stage Runge-Kutta algorithm, attaining second-order temporal precision for deterministic models. A sixth-order compact approach is utilized to tackle spatial discretization difficulties, providing great accuracy and processing efficiency. The strategy employs the Euler-Maruyama method to manage Wiener processes in stochastic models, facilitating flexibility to fluctuations in fluid flow dynamics. A numerical scheme has been proposed for solving stochastic time-dependent partial differential equations. The scheme was explicit, and it was comprised of two stages. The compact scheme was chosen to discretize space variable(s). The stability and consistency in the mean square sense of the scheme were presented. The scheme did not require any other scheme to get a solution to the problem. The proposed scheme is implemented in a dimensionless stochastic model of Williamson fluid dynamics, integrating essential physical phenomena such as Darcy–Forchheimer drag and the shear-thinning characteristics of the fluid. Numerical studies on the Stokes first issue indicate that the suggested technique surpasses current second-order Runge-Kutta methods in accuracy, especially in describing the complex relationship between deterministic and stochastic factors affecting the flow. The concluding points can be expressed as.

1. The modification of the proposed scheme performed better than the existing second-order Runge-Kutta method.

Finally, a computationally efficient and highly accurate method for non-Newtonian fluids in porous media is introduced, expanding the current state-of-the-art in stochastic fluid dynamics numerical techniques. This paradigm could be further developed in future studies to examine various types of non-Newtonian fluids and more complicated flow configurations, including three-dimensional geometries or transient boundary conditions. Furthermore, there is still room for improvement in the scheme’s optimization for real-time applications and large-scale simulations.