Numerical computational model for an unsteady hybrid nanofluid flow in a porous medium past an MHD rotating sheet

Salma Ahmedai; Precious Sibanda; Sicelo P. Goqo; Osman A.I. Noreldin

doi:10.1515/nleng-2022-0368

Article Open Access

Numerical computational model for an unsteady hybrid nanofluid flow in a porous medium past an MHD rotating sheet

Salma Ahmedai , Precious Sibanda , Sicelo P. Goqo and Osman A.I. Noreldin

Published/Copyright: September 2, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Nonlinear Engineering Volume 14 Issue 1

Abstract

This article presents a study of a rotating hybrid nanofluid flow in a porous medium past a stretching surface using an overlapping grid multi-domain bivariate spectral simple iteration method (OMD-BSSIM). The objective is to appraise the performance of the OMD-BSSIM on a system of nonlinear partial differential equations that describe the unsteady flow of a hybrid nanofluid composed of copper and aluminum oxide nanoparticles in water. Additionally, the study aims to evaluate the impact of different fluid and surface parameters on the behavior of the hybrid nanofluid. The improved OMD-BSSIM algorithm is implemented in MATLAB and tested for its convergence and accuracy. The results are compared with previously published findings that were obtained using the overlapping spectral local linearization method and the bivariate spectral relaxation method. The study explores the use of R line graphs as a visualization tool to identify patterns and relationships among the skin friction coefficients, local Nusselt number, and local Sherwood number of the hybrid nanofluid and parameters over time. The findings reveal several key trends and patterns, such as the positive correlation between the shrinking ratio parameter and the skin friction coefficient and the positive relationship between increasing the rotation parameter and an inverse relationship between an increase in the alumina volume fraction coefficient and the local Nusselt number of the hybrid nanofluid. This article provides valuable insights into the flow behavior of hybrid nanofluids and highlights the potential of alumina nanoparticles as a tool for enhancing the thermal and mechanical properties of the fluid.

Keywords: hybrid nanofluid; multi-domain overlapping grid; bivariate spectral simple iteration method

Nomenclature

x 1 , x 2 , x 3: Cartesian coordinate system
u 1 , u 2 , u 3: velocities in direction x 1 , x 2 , x 3
t: time
θ: hybrid nanofluid temperature
θ ∞: ambient temperature
θ 0: wall temperature
χ: hybrid nanofluid concentration
χ ∞: ambient concentration
χ 0: wall concentration
B 0: magnetic field coefficient
Q 0: heat generation/absorption coefficient
c p: specific heat capacity
g: acceleration due to gravity
a , b: shrinking parameters
f , h: velocity components after transformation
: Greek symbols
ρ: density
Γ: heat generation/absorption parameter
λ: porosity parameter
μ: dynamic viscosity
σ: electrical conductivity
Ω 0: angular velocity
Ω: rotation parameter
κ: thermal conductivity
γ: chemical reaction parameter
ϕ: nanoparticle volume fraction
δ: stretching/shrinking ratio parameter
η: coordinate for transformation space
ζ: coordinate for transformation time
ϑ: non-dimensional temperature
ψ: non-dimensional concentration
: Subscripts
s 1: alumina nanoparticle
s 2: copper nanoparticle
b f: base fluid
n f: nanofluid
h n f: hybrid nanofluid
: Abbreviations
S: total number of the grid points
T: vector transpose
BCs: boundary conditions
MHD: magnetohydrodynamic
SIM: simple iteration method
LLM: local linearization method
OMD: overlapping grid multi-domain
SSIM: spectral simple iteration method
SLLM: spectral local linearization method
BSRM: bivariate spectral relaxation method

1 Introduction

Hybrid nanofluids are a popular research topic due to, among other uses, their ability to improve heat transfer, thermal conductivity, and filtration efficiency [1]. Nanoparticles can be made from a variety of materials [2,3]. Adding nanoparticles to a liquid has been shown to improve electrical, thermal, and acoustical conductivity. Hybrid nanofluids consist of a base fluid and a mixture of nanoparticles [4]. The most commonly used base fluids include water, organic liquids, oils, bio-fluids, and polymeric solutions [5]. Numerous experiments and numerical simulations have demonstrated that nanofluids and hybrid nanofluids have high heat transfer efficiency compared to base fluids. However, researchers continue to push the boundaries of research in this area to expand knowledge. For example, some studies have explored the use of two different types of nanoparticles in nanofluids, rather than one type of particle. In addition, researchers have examined the impact of external factors, such as electric and magnetic fields, the effects of porous media, shrinking or stretching plates, and temperature, on the motion and heat transfer of hybrid nanofluids.

The field of magnetohydrodynamics (MHD) was first formally introduced by Alfvén [6]. MHD studies primarily focus on the dynamics and behavior of electrically conductive fluid, with particular emphasis on heat transfer management. The pioneering investigation of convective flow past a stretching sheet was conducted by Crane [7]. Many other studies have subsequently explored different aspects of the stretching sheet problem [8,9]. In manufacturing processes that involve stretching surfaces, it is crucial to take into account how the sheet is stretched and how heat within it is controlled. The dynamics of fluid motion through porous media has garnered significant attention. Understanding the nuances of fluid flow through porous media has ramifications for fields such as hydro-geology and petroleum engineering [10]. This requires an in-depth understanding of the intricate interplay among the fluid, the porous material, and the array of forces that govern and influence the flow through these intricate pathways. In this work, we combine the study of MHD effects with the impact of nanoparticles on a stretching surface. Despite the progress made in this fields, there is still a lot to learn about the mechanisms behind the enhanced heat transfer rates of hybrid nanofluids. This article aims to contribute to this understanding by studying the flow of hybrid nanofluids and examining the influence of external factors on their heat transfer properties.

Khalili et al. [11] performed a numerical investigation of the flow and heat transfer characteristics of nanofluids consisting of copper (Cu), titanium dioxide ( TiO 2 ), and aluminum oxide ( Al 2 O 3 ) in water ( H 2 O ) within a porous medium. Their study focused on stagnation point flow near stretching and shrinking sheets. They used a similarity transformation to reduce the governing nonlinear partial differential equations (N-PDEs) to a nonlinear ordinary differential equations (N-ODEs). They found that the permeability parameter had a more significant impact on the flow and heat transfer of a nanofluid than the magnetic parameter. Further, they observed that the skin friction and heat transfer coefficients increased with an increase in the nanoparticle volume fraction parameter. Dinarvand et al. [12] conducted a study of the unsteady two-dimensional flow of a hybrid nanofluid consisting of TiO 2 -Cu/H 2 O in a porous medium at the stagnation point of stretching and shrinking sheets. They employed the bvp4c numerical scheme, a MATLAB-based finite difference code, which uses the three-stage Lobatto IIIa formula [13] to solve the N-ODEs. They showed that the thermal characteristics of the hybrid nanofluid were higher than those of the base fluid (water only) and other single-component nanofluids. Their study revealed that as the volume fraction of Cu nanoparticles increased, the skin friction coefficient and local Nusselt number increased linearly. They noted that the heat transfer rate improved with increasing porosity and magnetic parameters, while it decreased with increasing unsteadiness parameter.

Khan et al. [14] recently investigated the behavior of ferrous oxide ( Fe 3 O 4 ), single-walled carbon nanotubes (SWCNTs), and multi-walled carbon nanotubes (MWCNTs) in water hybrid nanofluids between two parallel plates under a variable magnetic field. They employed two numerical schemes, the bvp4c and the parametric continuation method [15], to solve the transformed N-ODEs. They noted that an increase in the volume fraction of the nanoparticles led to a significant fluctuation in the velocity profile near the channel wall due to an increase in the fluid density. They observed that SWCNTs had a greater impact on temperature than MWCNTs. Mkhatshwa and Khumalo [16] examined the MHD two-dimensional unsteady flow over a stretching surface within a Darcy–Forchheimer porous medium. Their study focused on the flow of a temperature-variant hybrid nanofluid consisting of aluminum oxide and silver in water. They proposed a new overlapping grid spectral local linearization method (O-SLLM) [17] to solve the transformed N-ODEs. The authors found that fluid injection led to an increase in velocity, energy, and mass profiles. They observed that the rate of entropy generation could be minimized by reducing the porosity parameter and the Brinkman number, incorporating velocity slip conditions in the flow system, and using a shear-thinning Carreau fluid. This study offers valuable insights into the use of the overlapping grid approach and the factors that enable accurate results with a minimal number of grid points.

Despite the extensive investigations on hybrid nanofluids, to the best of our knowledge, there has been no attempt to find a numerical solution for the three-dimensional unsteady hybrid nanofluid flow. To address this gap, we model the MHD rotating flow of an Al 2 O 3 –Cu/H 2 O hybrid nanofluid in a porous medium past a shrinking surface with thermal effects. The transformed mathematical model is governed by a system of N-PDEs. This study aims to provide a comprehensive understanding of the behavior of hybrid nanofluids in three-dimensional unsteady flows and to highlight the factors that influence their flow and thermal properties. The numerical method employed in this study offers a promising approach for investigating the behavior of hybrid nanofluids in complex geometries and under different conditions. Researchers have employed various numerical methods to obtain solutions for N-PDEs that govern flow problems.

Recently, Mkhatshwa et al. [18,19] proposed an improvement of the standard algorithm for solving N-PDEs by introducing a novel overlapping grid multi-domain bivariate spectral quasi-linearization method (OMD-BSQLM) and the overlapping multi-domain bivariate spectral local linearization method (OMD-BSLLM). The solution algorithm involves linearizing and decoupling the N-PDEs using a uni-variate linearization technique along with a spectral collocation discretization. The methods have been successfully applied to solve various flow problems, including models for hybrid nanofluid flows. The combination of the overlapping technique with the bivariate spectral method has been demonstrated to provide accurate solutions [20]. The overlapping grid multi-domain method is used to solve N-PDEs involving two independent variables, namely space and time. The approach involves dividing the time domain into distinct, non-overlapping sub-domains, while the space domain is divided into sub-domains that overlap with each other. The overlapping techniques that have been developed guarantee a matrix with fewer elements and require a smaller number of grid points [17].

This article presents the overlapping grid multi-domain bivariate spectral simple iteration method (OMD-BSSIM) for finding numerical solutions to a system of N-PDEs that model unsteady hybrid nanofluid flow problems. The OMD-BSSIM is extended to the OMD-BSLLM by applying the simple iteration technique (SIM) instead of the local linearization technique (LLM). The method is modified to achieve computational efficiency. Linearization based on truncated Taylor series approximations is employed to simplify terms in the nonlinear differential equations [21]. The LLM, quasi-linearization method, and Keller-box methods are based on a one-term Taylor series expansion and are therefore susceptible to truncation errors. The relaxation method (RM) technique is based on the assumption that the non-linear terms are known from previous iterations. The SIM uses ideas akin to those of fixed-point iteration as an iterative scheme. Both the SIM and LLM techniques were proposed by Motsa et al. [22,23]. The effectiveness of the SIM is validated by comparison with the finite difference method for boundary-value problems [24].

2 Mathematical model

Consider the unsteady MHD flow of an incompressible, electrically conducting, rotating, stratified Cu- Al 2 O 3 ∕ H 2 O hybrid nanofluid flow past a linearly shrinking sheet. The coordinate system and the physical model of the flow problem are shown in Figure 1, where x 1 , x 2 , and x 3 represent the Cartesian coordinate system. The rotation of the hybrid nanofluid is along the vertical axis so that the fluid angular velocity Ω 0 is constant. The wall temperature θ 0 and the ambient temperature θ ∞ are positive constants. Under these assumptions, the continuity, momentum, energy, and concentration equations are given as follows [25,26]:

(1) ∂ u 1 ∂ x 1 + ∂ u 2 ∂ x 2 + ∂ u 3 ∂ x 3 = 0 ,

(2) ∂ u 1 ∂ t + u 1 ∂ u 1 ∂ x 1 + u 2 ∂ u 1 ∂ x 2 + u 3 ∂ u 1 ∂ x 3 − 2 Ω 0 u 2 = μ h n f ρ h n f − ∂ p ∂ x 1 + ∂ 2 u 1 ∂ x 3 2 − u 1 k p − σ h n f ρ h n f B 0 2 u 1 + g β h n f ( θ − θ ∞ ) ,

(3) ∂ u 2 ∂ t + u 1 ∂ u 2 ∂ x 1 + u 2 ∂ u 2 ∂ x 2 + u 3 ∂ u 2 ∂ x 3 + 2 Ω 0 u 1 = μ h n f ρ h n f − ∂ p ∂ x 2 + ∂ 2 u 2 ∂ x 3 2 − u 2 k p − σ h n f ρ h n f B 0 2 u 2 ,

(4) ∂ u 3 ∂ t + u 1 ∂ u 3 ∂ x 1 + u 2 ∂ u 3 ∂ x 2 + u 3 ∂ u 3 ∂ x 3 = μ h n f ρ h n f − ∂ p ∂ x 3 + ∂ 2 u 3 ∂ x 3 2 ,

(5) ∂ θ ∂ t + u 1 ∂ θ ∂ x 1 + u 2 ∂ θ ∂ x 2 + u 3 ∂ θ ∂ x 3 = α h n f ∂ 2 θ ∂ x 3 2 + Q w ( ρ c p ) h n f ( θ − θ ∞ ) ,

(6) ∂ χ ∂ t + u 1 ∂ χ ∂ x 1 + u 2 ∂ χ ∂ x 2 + u 3 ∂ χ ∂ x 3 = D h n f ∂ 2 χ ∂ x 3 2 + k r ( χ − χ ∞ ) ,

subject to the boundary conditions (BCs)

(7) t = 0 : u 1 = 0 , u 2 = 0 , u 3 = 0 , θ = θ 0 , χ = χ 0 ,

(8) t ≥ 0 : u 1 = U 1 , u 2 = U 2 , u 3 = 0 , θ = θ 0 , χ = χ 0 as x 3 = 0 ,

(9) t ≥ 0 : u 1 → 0 , u 2 → 0 , u 3 → 0 , θ → θ ∞ , χ → χ ∞ , as x 3 → ∞ .

The linearly shrinking velocities at the surface are defined by

U 1 = a x 1 , and U 2 = b x 2 ,

respectively. In Eqs. (1)–(6), μ h n f , ρ h n f , κ h n f , ( ρ c p ) h n f , σ h n f , and D h n f represent the dynamic viscosity, density, thermal conductivity, heat capacity, electrical conductivity, and diffusivity of the hybrid nanofluid. k p , Q 0 , B 0 , and k r are the permeability, heat generation or absorption, magnetic, and mass diffusion coefficients, respectively.

Figure 1

Geometry configuration of the problem.

Table 1 shows the thermophysical properties of H 2 O , Cu, and Al 2 O 3 at the reference temperature of 298 K. These values are given by Ahmed et al. [27] and Lund et al. [28]. The data in Table 2 have been experimentally validated in [29–31], where ϕ represents the volume fraction of nanoparticles.

Table 1

Alumina, copper, and water thermophysical properties [27,28]

Physical properties	H 2 O	Cu	Al 2 O 3
β ( K − 1 )	21 × 1 0 − 5	1.67 × 1 0 − 5	0.85 × 1 0 − 5
σ ( S m − 1 )	0.05	5.96 × 1 0 7	3.69 × 1 0 7
ρ ( kg m − 3 )	997.1	8,933	3,970
κ ( W m − 1 K − 1 )	0.613	400	40
c p ( J kg − 1 K − 1 )	4179	385	765

Table 2

Thermophysical features of nanofluid and hybrid nanofluid [29–31]

Properties	Nanofluid	Hybrid nanofluid
Density	( 1 − ϕ s 1 ) ρ b f + ϕ s 1 ρ s 1	( 1 − ϕ s 2 ) ρ n f + ϕ s 2 ρ s 2
Dynamic viscosity	μ b f ( 1 − ϕ s 1 ) 2.5	μ n f ( 1 − ϕ s 2 ) 2.5
Electrical conductivity	σ s 1 − 2 ϕ s 1 ( σ b f − σ s 1 ) + 2 σ b f σ s 1 + ϕ s 1 ( σ b f − σ s 1 ) + 2 σ b f σ b f	σ s 2 − 2 ϕ s 2 ( σ n f − σ s 2 ) + 2 σ n f σ s 2 + ϕ s 2 ( σ n f − σ s 2 ) + 2 σ n f σ n f
Heat capacity	( 1 − ϕ s 1 ) ( ρ c p ) b f + ϕ s 1 ( ρ c p ) s 1	( 1 − ϕ s 2 ) ( ρ c p ) n f + ϕ s 2 ( ρ c p ) s 2
Mass diffusivity	D b f ( 1 − ϕ s 1 ) 2.5	D n f ( 1 − ϕ s 2 ) 2.5
Thermal conductivity	κ s 1 − 2 ϕ s 1 ( κ b f − κ s 1 ) + 2 κ b f κ s 1 + ϕ s 1 ( κ b f − κ s 1 ) + 2 κ b f κ b f	κ s 2 − 2 ϕ s 2 ( κ n f − κ s 2 ) + 2 κ n f κ s 2 + ϕ s 2 ( κ n f − κ s 2 ) + 2 κ n f κ n f
Thermal expansion	( 1 − ϕ s 1 ) ( ρ β ) b f + ϕ s 1 ( ρ β ) s 1	( 1 − ϕ s 2 ) ( ρ β ) n f + ϕ s 2 ( ρ β ) s 2

Eqs. (1)–(6) are transformed into a non-dimensional form [25,32,33] using

(10) u 1 = a x 1 ∂ f ( ζ , η ) ∂ η , u 2 = a x 2 ∂ h ( ζ , η ) ∂ η , u 3 = − a ν b f ζ ( f ( ζ , η ) + h ( ζ , η ) ) , η = a ζ ν b f x 3 , ζ = 1 − e − τ , θ = θ ∞ + ( θ 0 − θ ∞ ) ϑ ( ζ , η ) , χ = χ ∞ + ( χ 0 − χ ∞ ) ψ ( ζ , η ) , τ = a t .

Eq. (1) is satisfied identically using the above similarity transformations. Substituting Eq. (10) into Eqs. (2)–(6) yields

(11) A 1 ∂ 3 f ∂ η 3 + ζ f + h − η 2 + η 2 ∂ 2 f ∂ η 2 − ζ ∂ f ∂ η + A 1 λ + A 2 M 2 ∂ f ∂ η + 2 Ω ζ ∂ h ∂ η + A 3 ζ ϑ = ( ζ − ζ 2 ) ∂ 2 f ∂ ζ ∂ η ,

(12) A 1 ∂ 3 h ∂ η 3 + ζ f + h − η 2 + η 2 ∂ 2 h ∂ η 2 − ζ ∂ h ∂ η + A 1 λ + A 2 M 2 ∂ h ∂ η − 2 Ω ζ ∂ f ∂ η = ( ζ − ζ 2 ) ∂ 2 h ∂ ζ ∂ η ,

(13) A 4 ∂ 2 ϑ ∂ η 2 + ζ f + h − η 2 + η 2 ∂ ϑ ∂ η + Γ ζ ϑ = ( ζ − ζ 2 ) ∂ ϑ ∂ ζ ,

(14) A 5 ∂ 2 ψ ∂ η 2 + ζ f + h − η 2 + η 2 ∂ ψ ∂ η − γ ζ ψ = ( ζ − ζ 2 ) ∂ ψ ∂ ζ .

Here, A 1 , A 2 , A 3 , A 4 , and A 5 are constants given by

A 1 = μ h n f ∕ μ b f ρ h n f ∕ ρ b f , A 2 = σ h n f ∕ σ b f ρ h n f ∕ ρ b f , A 3 = R i β h n f β b f , A 4 = 1 P r κ h n f ∕ κ b f ( ρ c p ) h n f ∕ ( ρ c p ) b f , A 5 = 1 S c D h n f D b f .

The BC Eqs. (7)–(9) became

(15) f ( ζ , 0 ) = 0 , h ( ζ , 0 ) = 0 , ∂ f ∂ η ( ζ , 0 ) = 1 , ∂ h ∂ η ( ζ , 0 ) = δ , ϑ ( ζ , 0 ) = 1 , ψ ( ζ , 0 ) = 1 , ζ ∈ [ 0 , 1 ] ,

(16) ∂ f ∂ η ( ζ , η ) → 0 , ∂ h ∂ η ( ζ , η ) → 0 , ϑ ( ζ , η ) → 0 , ψ ( ζ , η ) → 0 , ζ ∈ [ 0 , 1 ] , as η → ∞ ,

where δ is the stretching ( δ > 0 ) or shrinking ( δ < 0 ) parameter, when δ = 0 , the problem reduces to the two-dimensional case, λ is the porosity parameter, Ω is rotation parameter, M is the magnetic parameter, Γ is the heat absorption ( Γ < 0 ) and generation ( Γ > 0 ) parameter, and γ is the chemical reaction parameter, which are defined as

δ = b a , λ = ν b f a k p , Ω = Ω 0 a , M = B 0 σ b f a ρ b f , Γ = Q 0 a ( ρ c p ) b f , γ = k r a .

The parameter Ri is the local Richardson number, Gr is the local Grashof number, Re is the local Reynolds number, Pr is the Prandtl number, and Sc is the Schmidt number, which are defined as

Ri = Gr Re 2 , Gr = g Δ T β b f x 1 3 ν b f 2 , Re = U 0 x 1 ν b f , Pr = ν b f ( ρ c p ) b f κ b f , Sc = ν b f D b f .

The physical quantities of interest are skin friction coefficients along the x 1 and x 2 axis are C f x 1 , C f x 2 . We denote the local Nusselt number by Nu and the local Sherwood number by Sh, respectively [25]. Using (10), we obtain the dimensionless version of these quantities as

(17) Re x 1 1 2 ζ 1 2 C f x 1 = μ h n f μ b f ∂ 2 f ∂ η 2 ( ζ , 0 ) , Re x 2 1 2 ζ 1 2 C f x 1 = μ h n f μ b f ∂ 2 h ∂ η 2 ( ζ , 0 ) , Re x 1 − 1 2 ζ 1 2 Nu = − κ h n f κ b f ∂ ϑ ∂ η ( ζ , 0 ) , Re x 1 − 1 2 ζ 1 2 Sh = − D h n f D b f ∂ ψ ∂ η ( ζ , 0 ) .

3 Method of solution

In this section, we develop the iterative OMD-BSSIM scheme for the solution of the flow Eqs. (11)–(14) together with the BCs (15)–(16). The solution algorithm encompasses two fundamental components: the linearization of the N-PDEs utilizing the SIM technique and the incorporation of Chebyshev bivariate spectral collocation discretization.

3.1 SIM

In the SIM implementation, we consider the nonlinear terms as a combination of known (at iteration n ) and unknown (at iteration n + 1 ) functions during a specific iteration n . Consequently, in each nonlinear term of a given equation, the function with a higher derivative is treated as the unknown variable. By applying the SIM approach to the N-PDEs (11)–(14), we obtain a linear PDEs of the following form:

(18) α 10 , n ∂ 3 f n + 1 ∂ η 3 + α 11 , n ∂ 2 f n + 1 ∂ η 2 + α 12 , n ∂ f n + 1 ∂ η + α 13 , n f n + 1 + α 14 , n ∂ f n + 1 ∂ ζ + α 15 , n ∂ 2 f n + 1 ∂ ζ ∂ η = R 1 , n ,

(19) α 20 , n ∂ 3 h n + 1 ∂ η 3 + α 21 , n ∂ 2 h n + 1 ∂ η 2 + α 22 , n ∂ h n + 1 ∂ η + α 23 , n h n + 1 + α 24 , n ∂ h n + 1 ∂ ζ + α 25 , n ∂ 2 h n + 1 ∂ ζ ∂ η = R 2 , n ,

(20) α 30 , n ∂ 2 ϑ n + 1 ∂ η 2 + α 31 , n ∂ ϑ n + 1 ∂ η + α 32 , n ϑ n + 1 + α 33 , n ∂ ϑ n + 1 ∂ ζ = R 3 , n ,

(21) α 40 , n ∂ 2 ψ n + 1 ∂ η 2 + α 41 , n ∂ ψ n + 1 ∂ η + α 42 , n ψ n + 1 + α 43 , n ∂ ψ n + 1 ∂ ζ = R 4 , n .

In Eqs. (18)–(21), the coefficients α m r , n and the right-hand sides R m , n , where m = 1 , 2, 3, 4, are defined as

α 10 , n = A 1 , α 20 , n = A 1 , α 30 , n = A 4 , α 40 , n = A 5 , α 11 , n = ζ f n + h n − η 2 + η 2 , α 21 , n = ζ f n + 1 + h n − η 2 + η 2 , α 31 , n = ζ f n + 1 + h n + 1 − η 2 + η 2 , α 41 , n = ζ f n + 1 + h n + 1 − η 2 + η 2 , α 12 , n = − ζ ∂ f n ∂ η + A 1 λ + A 2 M 2 , α 22 , n = − ζ ∂ h n ∂ η + A 1 λ + A 2 M 2 , α 32 , n = Γ ζ , α 42 , n = − γ ζ , α 13 , n = 0 , α 23 , r = 0 , α 33 , n = ζ 2 − ζ , α 43 , n = ζ 2 − ζ , α 14 , n = 0 , α 24 , n = 0 , α 15 , n = ζ 2 − ζ , α 25 , n = ζ 2 − ζ , R 1 , n = − 2 Ω ∂ h n ∂ η + A 3 ϑ n , R 2 , n = 2 Ω ∂ f n + 1 ∂ η , R 3 , n = 0 , R 4 , n = 0 .

3.2 Chebyshev differentiation

The main feature of the OMD-Chebyshev spectral collocation method is the use of an overlapping grid procedure in the space variable ( η ). Accordingly, the space domain is decomposed into small overlapping subdomains of uniform length. On the contrary, the time variable ( ζ ) domain is partitioned into smaller non-overlapping sub-intervals [20]. Figure 2 shows the non-overlapping time domain I ∈ [ 0 , ζ F ] . The time domain is partitioned into p non-overlapping sub-intervals of equal length defined as

I ε = [ ζ ε − 1 , ζ ε ] , ζ ε − 1 < ζ ε , ε = 1 , 2 , 3 , 4 , … , p ,

$Figure 2 Non-overlapping grid ( ζ \zeta -domain).$

Figure 2

Non-overlapping grid ( ζ -domain).

Figure 3 displays the sub-division of the overlapping space domain ϒ ∈ [ 0 , η ∞ ] . The space domain is partitioned into q overlapping sub-intervals, which are specified by the following definition:

ϒ ϱ = [ η 0 ϱ , η N η ϱ ] , ϱ = 1 , 2 , 3 , … , q .

$Figure 3 Overlapping grid ( η \eta domain).$

Figure 3

Overlapping grid ( η domain).

We transform the sub-intervals I ε and ϒ ϱ into the interval [ − 1 , 1 ] through linear transformation mappings [19,20]. We assume that the solution can be approximated by a bivariate Lagrange interpolating polynomial of the form [34]

(22) f ( ε , ϱ ) ( ζ , η ) ≈ ∑ i = 0 N η ∑ j = 0 N ζ F ( ε , ϱ ) ( ζ j , η i ) L i ( η ) L j ( ζ ) , h ( ε , ϱ ) ( ζ , η ) ≈ ∑ i = 0 N η ∑ j = 0 N ζ H ( ε , ϱ ) ( ζ j , η i ) L i ( η ) L j ( ζ ) ,

(23) ϑ ( ε , ϱ ) ( ζ , η ) ≈ ∑ i = 0 N η ∑ j = 0 N ζ Θ ( ε , ϱ ) ( ζ j , η i ) L i ( η ) L j ( ζ ) , ψ ( ε , ϱ ) ( ζ , η ) ≈ ∑ i = 0 N η ∑ j = 0 N ζ Ψ ( ε , ϱ ) ( ζ j , η i ) L i ( η ) L j ( ζ ) .

The functions L i ( η ) and L j ( ζ ) are the characteristic Lagrange cardinal polynomials. From Eqs. (22) to (23), the first spatial derivatives are computed as

(24) ∂ f ( ε , ϱ ) ∂ η ( ζ ˆ i , η ˆ j ) = ∑ k = 0 N η D ̄ i , k F ( ε , ϱ ) ( ζ ˆ j , η ˆ k ) = D F j ( ε , ϱ ) , ∂ h ( ε , ϱ ) ∂ η ( ζ ˆ i , η ˆ j ) = ∑ k = 0 N η D ̄ i , k H ( ε , ϱ ) ( ζ ˆ j , η ˆ k ) = D H j ( ε , ϱ ) ,

(25) ∂ ϑ ( ε , ϱ ) ∂ η ( ζ ˆ i , η ˆ j ) = ∑ k = 0 N η D ̄ i , k Θ ( ε , ϱ ) ( ζ ˆ j , η ˆ k ) = D Θ j ( ε , ϱ ) , ∂ ψ ( ε , ϱ ) ∂ η ( ζ ˆ i , η ˆ j ) = ∑ k = 0 N η D ̄ i , k Ψ ( ε , ϱ ) ( ζ ˆ j , η ˆ k ) = D Ψ j ( ε , ϱ ) .

The matrix D i , k is the standard first-order Chebyshev differentiation matrix of size ( N η + 1 ) × ( N η + 1 ) . The D differentiation matrix represents the Chebyshev spatial differentiation matrix in the ϱ th sub-domain, and D is of size ( S + 1 ) × ( S + 1 ) , where S = N η + ( N η − 1 ) ( q − 1 ) is the total number of grid points in the whole η domain. The matrix D is defined as

(26) D = D 0,0 ( ε , q ) D 0 , 1 ( ε , q ) … D 0 , N η − 1 ( ε , q ) D 0 , N η ( ε , q ) D 1,0 ( ε , q ) D 1 , 1 ( ε , q ) … D 1 , N η − 1 ( ε , q ) D 1 , N η ( ε , q ) ⋱ ⋱ ⋱ ⋱ ⋱ D N η − 1,0 ( ε , q ) D N η − 1 , 1 ( ε , q ) … D N η − 1 , N η − 1 ( ε , q ) D N η − 1 , N η ( ε , q ) D 1,0 ( ε , q − 1 ) D 1 , 1 ( ε , q − 1 ) … D 1 , N η − 1 ( ε , q − 1 ) D 1 , N η ( ε , q − 1 ) D 2,0 ( ε , q − 1 ) D 2,1 ( ε , q − 1 ) … D 2 , N η − 1 ( ε , q − 1 ) D 2 , N η ( ε , q − 1 ) ⋱ ⋱ ⋱ ⋱ ⋱ D N η − 1,0 ( ε , q − 1 ) D N η − 1 , 1 ( ε , q − 1 ) … D N η − 1 , N η − 1 ( ε , q − 1 ) D N η − 1 , N η ( ε , q − 1 ) ⋱ ⋱ D 1,0 ( ε , 1 ) D 1 , 1 ( ε , 1 ) … D 1 , N η − 1 ( ε , 1 ) D 1 , N η ( ε , 1 ) D 2,0 ( ε , 1 ) D 2,1 ( ε , 1 ) … D 2 , N η − 1 ( ε , 1 ) D 2 , N η ( ε , 1 ) ⋱ ⋱ ⋱ ⋱ ⋱ D N η , 0 ( ε , 1 ) D N η , 1 ( ε , 1 ) … D N η , N η − 1 ( ε , 1 ) D N η , N η ( ε , 1 )

The time derivatives are computed as

(27) ∂ f ( ε , ϱ ) ∂ ζ ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι F ˆ j ( ε , ϱ ) , ∂ h ( ε , ϱ ) ∂ ζ ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι H ˆ j ( ε , ϱ ) , ∂ ϑ ( ε , ϱ ) ∂ ζ ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι Θ ˆ j ( ε , ϱ ) ,

(28) ∂ ψ ( ε , ϱ ) ∂ ζ ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι Ψ ˆ j ( ε , ϱ ) , ∂ 2 f ( ε , ϱ ) ∂ ζ ∂ η ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι D F ˆ j ( ε , ϱ ) , ∂ 2 h ( ε , ϱ ) ∂ ζ ∂ η ( ζ ˆ i , η ˆ j ) = ∑ ι = 0 N ζ d j , ι D H ˆ j ( ε , ϱ ) .

The matrix d ˆ j , ι is the standard first-order Chebyshev differentiation matrix of size ( N ζ + 1 ) × ( N ζ + 1 ) . The vector F ˆ j ( ε , ϱ ) is defined as

(29) F ˆ j ( ε , ϱ ) = [ f ( ε , ϱ ) ( ζ j ( ε ) , η 0 ( ϱ ) ) , f ( ε , ϱ ) ( ζ j ( ε ) , η 1 ( ϱ ) ) , f ( ε , ϱ ) ( ζ j ( ε ) , η 2 ( ϱ ) ) , … , f ( ε , ϱ ) ( ζ j ( ε ) , η N η ( ϱ ) ) ] T ,

where the superscript T representing vector transpose. Similarly for vectors H ˆ j ( ε , ϱ ) , Θ ˆ j ( ε , ϱ ) , and Ψ ˆ j ( ε , ϱ ) . The higher order derivatives with respect to the space variable obtained as

(30) ∂ s f ( ε , ϱ ) ∂ η s ( ζ ˆ i , η ˆ j ) = D s F j ( ε ) , ∂ s h ( ε , ϱ ) ∂ η s ( ζ ˆ i , η ˆ j ) = D s H j ( ε , ϱ ) , ∂ s ϑ ( ε , ϱ ) ∂ η s ( ζ ˆ i , η ˆ j ) = D s Θ j ( ε , ϱ ) , ∂ s ψ ( ε , ϱ ) ∂ η s ( ζ ˆ i , η ˆ j ) = D s Ψ j ( ε , ϱ ) .

Substituting Eqs. (22)–(30) into Eqs. (18)–(21), we obtain

(31) A ˆ F ˆ j , n + 1 = K ˆ 1 , j , n ,

(32) B ˆ H ˆ j , n + 1 = K ˆ 2 , j , n ,

(33) C ˆ Θ ˆ j , n + 1 = K ˆ 3 , j , n ,

(34) E ˆ Ψ ˆ j , n + 1 = K ˆ 4 , j , n .

Equations (31)–(34) represent the subsystem expressed as matrices of size N ζ ( S + 1 ) × N ζ ( S + 1 ) , where

A ˆ = α ˆ 10 , n ( ε , ϱ ) D 3 + α ˆ 11 , n ( ε , ϱ ) D 2 + α ˆ 12 , n ( ε , ϱ ) D + α ˆ 13 , n ( ε , ϱ ) I , K ˆ 1 , j , n = R ˆ 1 , j , n ( ε , ϱ ) − α ˆ 14 , n ( ε , ϱ ) d j , N ζ D F ˆ N ζ , n ( ε , ϱ ) , B ˆ = α ˆ 20 , n ( ε , ϱ ) D 3 + α ˆ 21 , n ( ε , ϱ ) D 2 + α ˆ 22 , n ( ε , ϱ ) D + α ˆ 23 , n ( ε , ϱ ) I , K ˆ 2 , j , n = R ˆ 2 , j , n ( ε , ϱ ) − α ˆ 24 , n ( ε , ϱ ) d j , N ζ D H ˆ N ζ , n ( ε , ϱ ) , C ˆ = α ˆ 30 , n ( ε , ϱ ) D 2 + α ˆ 31 , n ( ε , ϱ ) D + α ˆ 32 , n ( ε , ϱ ) I , K ˆ 3 , j , n = R ˆ 3 , j , n ( ε , ϱ ) − α ˆ 33 , n ( ε , ϱ ) d j , N ζ Θ ˆ N ζ , n ( ε , ϱ ) , E ˆ = α ˆ 40 , n ( ε , ϱ ) D 2 + α ˆ 41 , n ( ε , ϱ ) D + α ˆ 42 , n ( ε , ϱ ) I , K ˆ 4 , J , n = R ˆ 4 , j , n ( ε , ϱ ) − α ˆ 43 , n ( ε , ϱ ) d j , N ζ Ψ ˆ N ζ , n ( ε , ϱ ) ,

and I is the identity matrix of size N ζ ( S + 1 ) × N ζ ( S + 1 ) . The BC equations (15)–(16)

(35) F n + 1 ( ε , ϱ ) ( ζ j , η S ) = 0 , ∑ k = 0 S D S , k F n + 1 ( ε , ϱ ) ( ζ j , η k ) = 1 , ∑ k = 0 S D 0 , k F n + 1 ( ε , ϱ ) ( ζ j , η k ) = 0 ,

(36) H n + 1 ( ε , ϱ ) ( ζ j , η S ) = 0 , ∑ k = 0 S D S , k H n + 1 ( ε , ϱ ) ( ζ j , η k ) = δ , ∑ k = 0 S D 0 , k H n + 1 ( ε , ϱ ) ( ζ j , η k ) = 0 ,

(37) Θ n + 1 ( ε , ϱ ) ( ζ j , η S ) = 1 , Θ n + 1 ( ε , ϱ ) ( ζ j , η 0 ) = 0 , Ψ n + 1 ( ε , ϱ ) ( ζ j , η S ) = 1 , Ψ n + 1 ( ε , ϱ ) ( ζ j , η 0 ) = 0 .

The functions taken as guesses for initiating the iteration procedure are

(38) F ˆ j , 0 ( ε , ϱ ) = 1 − e − η i , H ˆ j , 0 ( ε , ϱ ) = δ ( 1 − e − η i ) , Θ ˆ j , 0 ( ε , ϱ ) = e − η i , Ψ ˆ j , 0 ( ε , ϱ ) = e − η i ,

which are chosen to satisfy the BC equations (35)–(37). After imposing the BC equations (35)–(37) into the matrix subsystems (31)–(34), the solutions are derived through solving the matrix systems iteratively

(39) F ˆ j , n + 1 = inv ( A ˆ ) K ˆ 1 , j , n ,

(40) H ˆ j , n + 1 = inv ( B ˆ ) K ˆ 2 , j , n ,

(41) Θ ˆ j , n + 1 = inv ( C ˆ ) K ˆ 3 , j , n ,

(42) Ψ ˆ j , n + 1 = inv ( E ˆ ) K ˆ 4 , j , n .

where inv computes the inverse of the matrix.

It can be seen that when τ → ∞ , the NL-PDEs (11)–(14) become

(43) A 1 f ‴ + ( f + h ) f ″ − ( A 1 λ + A 2 M 2 ) f ′ − f ′ 2 + 2 Ω h ′ + A 3 ϑ = 0 ,

(44) A 1 h ‴ + ( f + h ) h ″ − ( A 1 λ + A 2 M 2 ) h ′ − h ′ 2 − 2 Ω f ′ = 0 ,

(45) A 4 ϑ ″ + ( f + h ) ϑ ′ + Γ ϑ = 0 ,

(46) A 5 ψ ″ + ( f + h ) ψ ′ − γ ψ = 0 .

Applying the overlapping grid on the spatial variable ( η ) only, Eqs (43)–(46) can be solved using the O-SLLM [16] and the overlapping grid spectral simple iteration method (O-SSIM).

4 Results and discussion

To verify the accuracy of the SIM method, we utilized the O-SLLM and O-SSIM to solve Eqs. (43)–(46). Table 3 presents a comparison and validation demonstrating the clear consistency of the SIM. We used η ∞ = 10 , q = 12 , and N η = 12 . The accuracy of the numerical solution is confirmed through comparison with recently published results obtained using O-SLLM [16], as well as earlier published studies [35–37] for the case where R i = ϕ s 1 = ϕ s 2 = γ = M = Ω = Γ = 0 . It is noted that the O-SSIM gives accurate results for various Prandtl numbers in a few iterations n = 10 . Table 3 shows a comparison of the computational time (CPU time) for both methods to give converged results. It is noted that the O-SSIM approach converges faster than the O-SLLM approach.

Table 3

Comparison of overlapping SSIM for values of − ϑ ′ ( 0 ) at ζ = 1

Pr	Devi and Devi [35]	Khashii’ie et al. [36]	Waini et al. [37]	O-SLLM [16]	O-SSIM
2.00	0.91135	0.91135	0.911357	0.91135276	0.9113527637
6.13	1.75968	1.75968	1.759682	1.75968170	1.7596817080
7.00	1.89540	1.89540	1.895400	1.89540040	1.8954004012
20.0	3.35390	3.35390	3.353893	3.35390184	3.3538961832
CPU time (s)	—	—	—	0.052554	0.041551

Equations (11)–(14) were solved numerically using the OMD-BSSIM. Table 4 shows f ″ ( ζ , 0 ) , h ″ ( ζ , 0 ) , ϑ ′ ( ζ , 0 ) , and ψ ′ ( ζ , 0 ) for different values of ζ F using p = 20 , and N ζ = 10 in the time domain. To determine the accuracy of the OMD-BSSIM, numerical results are compared with Magagula et al. [33]. Table 4 gives a comparison of OMD-BSSIM for η ∞ = 20 , q = 5 , and N η = 12 . It is observed that the OMD-BSSIM gives better results than BSRM.

Table 4

Comparison of implementation the OMD-BSSIM when Ω = 0 , Ri = 0 , ϕ s 1 = 0 , ϕ s 2 = 0 , Sc = 1 , Pr = 1.5 , M = 2 , δ = 0.5 , λ = 0.5 , and γ = 1

	f ″ ( ζ , 0 )		h ″ ( ζ , 0 )		θ ′ ( ζ , 0 )		ψ ′ ( ζ , 0 )
ζ F	BSRM [33]	Present	BSRM [33]	Present	BSRM [33]	Present	BSRM [33]	Present
0.1	− 0.851257	− 0.851270	− 0.417150	− 0.417158	− 0.710882	− 0.710809	− 0.634443	− 0.634466
0.3	− 1.316705	− 1.316736	− 0.639602	− 0.639620	− 0.742842	− 0.742783	− 0.766867	− 0.766880
0.5	− 1.685306	− 1.685591	− 0.817649	− 0.817899	− 0.765244	− 0.765198	− 0.891207	− 0.891211
0.7	− 1.992608	− 1.993602	− 0.966603	− 0.967595	− 0.777270	− 0.777225	− 1.010045	− 1.010034
0.9	− 2.259335	− 2.261500	− 1.095983	− 1.098352	− 0.770807	− 0.770733	− 1.125549	− 1.125508

The numerical computations are executed using the values

(47) ϕ s 1 = 0.1 , ϕ s 2 = 0.02 , Ri = 3 , Pr = 6.9 , Sc = 2.22 , δ = − 1 , Ω = 1.5 , M 2 = 4 ,

(48) γ = 1.5 , Γ = 1.5 , λ = 2 , n = 40 , p = 20 N ζ = 5 , η ∞ = 10 , ζ F = 0.94 ,

unless indicated otherwise. We examined the residual errors and the solution errors for different values of q , N η , and S . We define the operators Λ m where m = 1 , 2, 3, 4, to represent the N-PDEs (11)–(14). The approximated solutions are denoted by F , H , Θ , and Ψ . The residual errors are defined as

Res f = ∣ Λ 1 ( F j , n , H j , n , Θ j , n , Ψ j , n ) ∣ , Res h = ∣ Λ 2 ( F j , n + 1 , H j , n , Θ j , n , Ψ j , n ) ∣ , Res ϑ = ∣ Λ 3 ( F j , n + 1 , H j , n + 1 , Θ j , n , Ψ j , n ) ∣ , Res ψ = ∣ Λ 4 ( F j , n + 1 , H j , n + 1 , Θ j , n + 1 , Ψ j , n ) ∣ .

Table 5 and Figure 4 present the residual errors for different values of q , N η , and S . We note in Table 5 that the OMD-BSSIM with q = 4 and N η = 15 gives the most accurate results and requires the least computational time. We note that the accuracy further improves by increasing q . Figure 4 indicates the variation of the residual errors for q = 4 and N η = 15 in both the spatial and temporal domains. Figures 4(a) and (b) illustrate the maximum residual errors Res f and Res h , respectively, revealing comparable magnitudes at approximately O ( 1 0 − 13 ) . In contrast, Figures 4(c) and (d) show the maximum residual errors of approximately O ( 1 0 − 14 ) for Res ϑ and Res ψ . The incorporation of solution updating strategies has significantly contributed to the improved accuracy. We defined the solution errors as

E f = ∥ F j , n − F j , n − 1 ∣ ∞ , E h = ∥ H j , n − H j , n − 1 ∣ ∞ , E ϑ = ∥ Θ j , n − Θ j , n − 1 ∣ ∞ , E ψ = ∥ Ψ j , n − Ψ j , n − 1 ∣ ∞ .

Table 6 illustrates how increasing the number of overlapping sub-intervals impacts the errors in the solution. We remark that when N η is less than q , the method gives better solutions errors. The overlapping technique generates accurate results using a relatively small N η throughout the spatial domain.

Table 5

Residual errors at ζ F = 0.94

N η	q	S	∥ Res f ∣ ∞	∥ Res h ∣ ∞	∥ Res ϑ ∣ ∞	∥ Res ψ ∣ ∞	CPU time in seconds
120	1	120	1.24792449 × 1 0 − 7	1.01281581 × 1 0 − 7	2.14038232 × 1 0 − 11	7.40120673 × 1 0 − 10	1.004388
20	5	96	4.38178501 × 1 0 − 10	4.84201141 × 1 0 − 10	2.21156427 × 1 0 − 13	1.78995707 × 1 0 − 12	0.891982
9	9	73	1.21587863 × 1 0 − 11	1.08099463 × 1 0 − 11	5.29576383 × 1 0 − 14	3.36512068 × 1 0 − 13	0.829143
4	15	46	8.38400584 × 1 0 − 13	2.48603492 × 1 0 − 13	2.22044605 × 1 0 − 14	2.55212518 × 1 0 − 14	0.767800

$Figure 4 The residual error graphs of the OMD-BSSIM solution when q = 15 q=15 , and N η = 4 {N}_{\eta }=4 . (a) Residual error graph of Λ 1 {\Lambda }_{1} , (b) residual error graph of Λ 2 {\Lambda }_{2} , (c) residual error graph of Λ 3 {\Lambda }_{3} , and (d) residual error graph of Λ 4 {\Lambda }_{4} .$

Figure 4

The residual error graphs of the OMD-BSSIM solution when q = 15 , and N η = 4 . (a) Residual error graph of Λ 1 , (b) residual error graph of Λ 2 , (c) residual error graph of Λ 3 , and (d) residual error graph of Λ 4 .

Table 6

Convergence of the solutions F , H , Θ , and Ψ at the iteration n = 40 , and ζ F = 0.94

N η	q	S	E f	E h	E ϑ	E ψ
120	1	120	9.17257781 × 1 0 − 9	1.37713215 × 1 0 − 8	1.86462135 × 1 0 − 9	5.15595289 × 1 0 − 10
20	5	96	2.09812168 × 1 0 − 10	3.25510729 × 1 0 − 11	8.16990919 × 1 0 − 12	1.96342942 × 1 0 − 12
9	9	73	2.24931185 × 1 0 − 12	1.58584257 × 1 0 − 12	2.36699549 × 1 0 − 13	1.12909682 × 1 0 − 13
4	15	46	1.12354570 × 1 0 − 13	7.86037901 × 1 0 − 14	7.37188088 × 1 0 − 14	4.14529522 × 1 0 − 14

Figure 5 illustrates the numerical values of velocity, temperature, and concentration functions when N η = 4 , and q = 15 .

$Figure 5 The numerical solution utilized the OMD-BSSIM. (a) Numerical solution of f ′ ( ζ , η ) f^{\prime} \left(\zeta ,\eta ) , (b) numerical solution of h ′ ( ζ , η ) h^{\prime} \left(\zeta ,\eta ) , (c) numerical solution of ϑ ( ζ , η ) {\vartheta }\left(\zeta ,\eta ) , and (d) numerical solution of ϕ ( ζ , η ) \phi \left(\zeta ,\eta ) .$

Figure 5

The numerical solution utilized the OMD-BSSIM. (a) Numerical solution of f ′ ( ζ , η ) , (b) numerical solution of h ′ ( ζ , η ) , (c) numerical solution of ϑ ( ζ , η ) , and (d) numerical solution of ϕ ( ζ , η ) .

We investigated the effect of the shrinking and stretching parameter ( δ ), the porosity parameter ( λ ), the magnetic parameter ( M ), the volume fraction of aluminum oxide nanoparticles ( ϕ s 1 ), the volume fraction of copper nanoparticles ( ϕ s 2 ), the rotation parameter ( Ω ), the chemical reaction parameter ( γ ), and the Schmidt number (Sc) on the fluid velocities, temperature, and concentration profiles of the unsteady Al 2 O 3 -Cu ∕ H 2 O hybrid nanofluid flow in porous medium when q = 8 , and N η = 20 .

Figure 6 shows the influence of the shrinking and stretching parameter, the porosity parameter, the magnetic parameter, and the volume fraction of aluminum oxide nanoparticle parameters on the velocity in the x 1 direction. Figure 6(a) shows that velocity tends to initially increase before decreasing as the shrinking parameter decreases, but it decreases in the case of the stretching parameter. Reducing the shrinking parameter ( a < b ) reduces the deformation experienced by the nanoparticles, resulting in an increase in the effective viscosity of the hybrid nanofluid, which causes an initial increase in velocity followed by a decrease further downstream in the x 1 direction. On the other hand, increasing the stretching parameter ( b > a ) leads to higher deformation of the fluid elements, which causes higher shear stress and an immediate decrease in velocity. Figure 6(b) indicates that increasing the porosity parameter improves the velocity by enhancing the permeability of the porous medium. A higher porous space means that there are more open channels for the fluid to flow through, resulting in less drag and a higher velocity. Figure 6(c) shows that an increase in the magnetic parameter value leads to a decline in the flow due to the resistive force, also known as the Lorentz force, increasing. Figure 6(d) shows that velocity increases with an increase in the volume fraction of aluminum oxide nanoparticles. This is because an increase in the solid particle volume fraction in a nanofluid leads to an increase in its effective viscosity due to the presence of additional particles in the fluid. This increase in viscosity results in an increase in the drag force experienced by the fluid as it flows, leading to a corresponding increase in the velocity. The influence of the shrinking and stretching parameter, the rotation parameter, the volume fraction of aluminum oxide nanoparticles, and the volume fraction of copper nanoparticles on the velocity profile in the x 2 direction is shown in Figure 7. Figure 7(a) shows that velocity increases with the shrinking parameter but decreases with the stretching parameter. It is noted that shrinking leads to a dominance of reverse flow, similar to the velocity profile in the case shown in Figure 6(a). Figure 7(b) displays the effect of varying the rotation parameter. It is observed that an increase in the rotation parameter can lead to an increase in velocity in the x 2 direction. This is because rotation can induce additional fluid mixing and enhance the transport of momentum in the fluid. The prevailing reverse flow is also observed. Figure 7(c) shows the impact of increasing the volume fraction of aluminum oxide nanoparticles on the velocity. It is noted that the flow velocity increases significantly and dominates the reverse flow. Figure 7(d) shows that the velocity in the x 2 direction decreases when the volume fraction of copper nanoparticles increases. This is because an increase in copper nanoparticles concentration leads to an increase in viscosity, causing an increase in flow resistance.

$Figure 6 Velocity profiles in the x 1 {x}_{1} direction at ζ = 0.94 \zeta =0.94 . (a) Fluctuation of f ′ ( ζ , η ) f^{\prime} \left(\zeta ,\eta ) along with δ \delta , (b) fluctuation of f ′ ( ζ , η ) f^{\prime} \left(\zeta ,\eta ) along with λ \lambda , (c) fluctuation of f ′ ( ζ , η ) f^{\prime} \left(\zeta ,\eta ) along with M M , and (d) fluctuation of f ′ ( ζ , η ) f^{\prime} \left(\zeta ,\eta ) along with ϕ s 1 {\phi }_{s1} .$

Figure 6

Velocity profiles in the x 1 direction at ζ = 0.94 . (a) Fluctuation of f ′ ( ζ , η ) along with δ , (b) fluctuation of f ′ ( ζ , η ) along with λ , (c) fluctuation of f ′ ( ζ , η ) along with M , and (d) fluctuation of f ′ ( ζ , η ) along with ϕ s 1 .

$Figure 7 Velocity profiles in the x 2 {x}_{2} direction at ζ = 0.94 \zeta =0.94 . (a) Fluctuation of h ′ ( ζ , η ) h^{\prime} \left(\zeta ,\eta ) along with δ \delta , (b) fluctuation of h ′ ( ζ , η ) h^{\prime} \left(\zeta ,\eta ) along with Ω \Omega , (c) fluctuation of h ′ ( ζ , η ) h^{\prime} \left(\zeta ,\eta ) along with ϕ s 1 {\phi }_{s1} , and (d) fluctuation of h ′ ( ζ , η ) h^{\prime} \left(\zeta ,\eta ) along with ϕ s 2 {\phi }_{s2} .$

Figure 7

Velocity profiles in the x 2 direction at ζ = 0.94 . (a) Fluctuation of h ′ ( ζ , η ) along with δ , (b) fluctuation of h ′ ( ζ , η ) along with Ω , (c) fluctuation of h ′ ( ζ , η ) along with ϕ s 1 , and (d) fluctuation of h ′ ( ζ , η ) along with ϕ s 2 .

Figure 8 shows the effect of the shrinking and stretching parameter, the porosity parameter, the rotation parameter, and the volume fraction of aluminum oxide nanoparticles on temperature. Figure 8(a) depicts that the temperature of the hybrid nanofluid flow increases with the shrinking parameter but decreases with the stretching parameter. The temperature of the hybrid nanofluid flow increases due to heat conduction from the surface, which can cause changes in density and viscosity. Figure 8(b) indicates that increasing the porosity parameter improves temperature distribution. Higher porosity can lead to lower thermal conductivity and viscosity of the fluid, resulting in reduced heat transfer. However, it should be noted that the presence of nanoparticles can enhance heat transfer by increasing the effective thermal conductivity and convective heat transfer coefficient of the hybrid nanofluid. Nevertheless, if the nanoparticles agglomerate and form blockages, they can reduce the flow rate, leading to a higher fluid temperature. Figure 8(c) shows that the temperature profile increases with the introduction of the rotation parameter. This is because increased rotation can induce higher shear rates in the hybrid nanofluid, leading to increased viscous dissipation and heating of the fluid. Figure 8(d) shows that temperature improves as the value of the volume fraction of aluminum oxide nanoparticles increases. This increase in thermal conductivity can reduce the thermal boundary layer and improve the heat transfer rate.

$Figure 8 Temperature profiles at ζ = 0.94 \zeta =0.94 . (a) Fluctuation of ϑ ( ζ , η ) {\vartheta }\left(\zeta ,\eta ) along with δ \delta , (b) fluctuation of ϑ ( ζ , η ) {\vartheta }\left(\zeta ,\eta ) along with λ \lambda , (c) fluctuation of ϑ ( ζ , η ) {\vartheta }\left(\zeta ,\eta ) along with Ri, and (d) fluctuation of ϑ ( ζ , η ) {\vartheta }\left(\zeta ,\eta ) along with ϕ s 1 {\phi }_{s1} .$

Figure 8

Temperature profiles at ζ = 0.94 . (a) Fluctuation of ϑ ( ζ , η ) along with δ , (b) fluctuation of ϑ ( ζ , η ) along with λ , (c) fluctuation of ϑ ( ζ , η ) along with Ri, and (d) fluctuation of ϑ ( ζ , η ) along with ϕ s 1 .

The effect of the shrinking and stretching parameter, the chemical reaction parameter, the Schmidt number, and the volume fraction of copper nanoparticles on the concentration profile, is illustrated in Figure 9. Figure 9(a) shows that the concentration of the hybrid nanofluid increases with the shrinking of the boundary layer and decreases with the stretching of the boundary layer. The concentration profile becomes more uniform due to increased mixing caused by turbulence and the thickness of the boundary layer. Figure 9(b) demonstrates that increasing chemical reaction parameter leads to a reduction in concentration. This is because the reaction will consume the nanoparticles in regions where the reaction is more active, leading to a decrease in the local concentration of nanoparticles. This effect can be more significant in regions where the reaction rate is high. Figure 9(c) reveals that the concentration decreases with an increase in the Schmidt number value. The Schmidt number is a dimensionless parameter that represents the ratio of momentum diffusivity to mass diffusivity in a fluid. It describes the transport of a passive scalar, such as the concentration of nanoparticles in a hybrid nanofluid. Increasing the Schmidt number leads to decreased mixing and slower transport of nanoparticles in the fluid. Figure 9(d) indicates a rising trend in concentration with increasing values of volume fraction of copper nanoparticles. However, it can be seen that the change in concentration is small.

$Figure 9 Concentration profiles at ζ = 0.94 \zeta =0.94 . (a) Fluctuation of ψ ( ζ , η ) \psi \left(\zeta ,\eta ) along with δ \delta , (b) fluctuation of ψ ( ζ , η ) \psi \left(\zeta ,\eta ) along with γ \gamma , (c) fluctuation of ψ ( ζ , η ) \psi \left(\zeta ,\eta ) along with Sc, and (d) fluctuation of ψ ( ζ , η ) \psi \left(\zeta ,\eta ) along with ϕ s 2 {\phi }_{s2} .$

Figure 9

Concentration profiles at ζ = 0.94 . (a) Fluctuation of ψ ( ζ , η ) along with δ , (b) fluctuation of ψ ( ζ , η ) along with γ , (c) fluctuation of ψ ( ζ , η ) along with Sc, and (d) fluctuation of ψ ( ζ , η ) along with ϕ s 2 .

This study also aims to examine how the copper and aluminum oxide fraction coefficients, shrinking ratio, rotation, and magnetic parameters relate to the skin friction coefficients ( C f x 1 , C f x 2 ), the local Nusselt number (Nu) and the local Sherwood number (Sh). To achieve this objective, we employed R visualization methods, specifically utilizing line graphs to visualize the trends and patterns (Figure 10). The line graphs represent the relationships between several dependent variables ( C f x 1 , C f x 2 , Nu, and Sh) and a range of parameters ( M , δ , ϕ s 1 , and ϕ s 2 ) over time. The horizontal axis corresponds to time, which ranges from 0.047 to 0.94. The vertical axis represents the values of the dependent variables, which are calculated based on changes in the parameters. Compared to other parameters over time, the graph displays an inverse relationship between the magnetic parameter and the skin friction coefficient in the x 1 direction, and a linear relationship between the magnetic parameter and the skin friction coefficient in the x 2 direction, while the effect of increasing stretching parameter appears to be completely opposite, with higher values of the stretching parameter corresponding to higher values of the skin friction coefficient in the x 1 direction and lower values of the skin friction coefficient in the x 2 direction. Also, we noted a consistent linear relationship between the volume fraction of aluminum oxide nanoparticles and the local Sherwood number, as well as between the rotation parameter and the local Nusselt number, with higher values of these parameters associated with higher values of the respective dependent variables. The values of the magnetic parameter, the volume fraction of copper nanoparticles, and the stretching parameter, do not seem to have a strong or direct effect on the values of Sherwood number or Nusselt number, at least not compared to the effect of volume fraction of aluminum oxide nanoparticles or rotation parameter.

$Figure 10 OMD-SSIM for C f x 1 {C}_{fx1} , C f x 2 {C}_{fx2} , Nu, and Sh when Pr = 6.90 {\rm{\Pr }}=6.90 .$

Figure 10

OMD-SSIM for C f x 1 , C f x 2 , Nu, and Sh when Pr = 6.90 .

5 Conclusion

The study presented the development of the OMD-BSSIM to solve model equations for the unsteady three-dimensional MHD flow of an incompressible, electrically conducting, rotating Cu- Al 2 O 3 ∕ H 2 O hybrid nanofluid flow past a linearly shrinking sheet. The solution algorithm combines the linearization of N-PDEs using the SIM technique and the use of the Chebyshev bivariate spectral collocation discretization. The method was shown to be effective in accurately capturing flow characteristics and achieving convergence. The method developed in this study can be used to solve systems of highly nonlinear partial differential equations that arise from a variety of real-world phenomena, including engineering problems. The main findings of this study are as follows:

The OMD-BSSIM may be used in both small and large computational domains.
The SIM requires less computational time than the LLM.
The incorporation of the overlapping grid approach contributes to improving the efficiency of the Chebyshev bivariate spectral collocation discretization.
The accuracy improves by increasing the number of the overlapping sub-intervals.
The OMD-BSSIM gives better solutions errors using a minimal number of grid points.
The analysis of residual errors highlights the significance of solution-updating strategies in improving the accuracy and convergence of the method.

The study found that increasing the nanoparticle aluminum oxide volume fraction and the shrinking parameter enhances the velocity and temperature, and the incremental increase in porous space improves the temperature of the hybrid nanofluid. Additionally, the study found that the Nusselt number increases with the stratified parameter, and that the Schmidt number is reduced by increasing the magnetic parameter. These findings provide valuable insights into the dynamics of hybrid nanofluids and highlight the potential of alumina nanoparticles as a means of enhancing the thermal and mechanical properties of the fluid. The study explored the use of R line graphs as a visualization tool to identify patterns and relationships between the skin friction coefficients, local Nusselt number, and Sherwood number of the hybrid nanofluid and parameters over time. The results suggest an inverse relationship between the strength of the magnetic and skin friction coefficient, and a linear relationship between rotation and the Nusselt number.

sahmedai@neelain.edu.sd

Acknowledgments

The authors are grateful to the University of KwaZulu-Natal.

Funding information: Authors state no funding involved.
Author contributions: Salma Ahmedai: conceptualization, methodology, software, investigation, writing – original draft preparation. Precious Sibanda: validation, writing – review and editing, visualization, supervision. Sicelo P. Goqo: writing-review and editing, visualization, supervision. Osman A.I. Noreldin: writing-review and editing, visualization, supervision. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: Authors state no conflict of interest.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Xu HJ, Xing ZB, Wang F, Cheng Z. Review on heat conduction, heat convection, thermal radiation and phase change heat transfer of nanofluids in porous media: Fundamentals and applications. Chem Eng Sci. 2019;195:462–83. 10.1016/j.ces.2018.09.045Search in Google Scholar

[2] Sarkar J, Ghosh P, Adil A. A review on hybrid nanofluids: recent research, development and applications. Renew Sustain Energy Rev. 2015;43:164–77. 10.1016/j.rser.2014.11.023Search in Google Scholar

[3] Muneeshwaran M, Srinivasan G, Muthukumar P, Wang CC. Role of hybrid-nanofluid in heat transfer enhancement-A review. Int Commun Heat Mass Transfer. 2021;125:105341. 10.1016/j.icheatmasstransfer.2021.105341Search in Google Scholar

[4] Li Y, Zhou J, Tung S, Schneider E, Xi S. A review on development of nanofluid preparation and characterization. Powder Tech. 2009;196(2):89–101. 10.1016/j.powtec.2009.07.025Search in Google Scholar

[5] Jahar S. A critical review on convective heat transfer correlations of nanofluids. Renew Sustain Energy Reviews. 2011;15(6):3271–7. 10.1016/j.rser.2011.04.025Search in Google Scholar

[6] Alfvén H. Existence of electromagnetic-hydrodynamic waves. Nature. 1942;150(3805):405–6. 10.1038/150405d0Search in Google Scholar

[7] Crane LJ. Flow past a stretching plate. Zeitschrift für angewandte Mathematik und Physik ZAMP. 1970;21:645–7. 10.1007/BF01587695Search in Google Scholar

[8] Chen CH. Effect of viscous dissipation on heat transfer in a non-Newtonian liquid film over an unsteady stretching sheet. J Non-Newton Fluid Mech. 2006;135(2–3):128–35. 10.1016/j.jnnfm.2006.01.009Search in Google Scholar

[9] Awan AU, Aziz M, Ullah N, Nadeem S, Abro KA. Thermal analysis of oblique stagnation point flow with slippage on second-order fluid. J Thermal Anal Calorimetry. 2021:147:3839–51. 10.1007/s10973-021-10760-zSearch in Google Scholar

[10] Yarlagadda A, Yoganathan A. Experimental studies of model porous media fluid dynamics. Experiments Fluids. 1989;8(1–2):59–71. 10.1007/BF00203066Search in Google Scholar

[11] Khalili S, Dinarvand S, Hosseini R, Tamim H, Pop I. flow and heat transfer near stagnation point over a stretching/shrinking sheet in porous medium filled with a nanofluid. Chin Phys B. 2014;23(4):048203. 10.1088/1674-1056/23/4/048203Search in Google Scholar

[12] Dinarvand S, Yousefi M, Chamkha A, et al. Numerical simulation of unsteady flow toward a stretching/shrinking sheet in porous medium filled with a hybrid nanofluid. J Appl Comput Mech. 2019;8:11–20. Search in Google Scholar

[13] Kierzenka J, Shampine LF. A BVP solver based on residual control and the Maltab PSE. ACM Trans Math Software (TOMS). 2001;27(3):299–316. 10.1145/502800.502801Search in Google Scholar

[14] Khan MS, Mei S, Shabnam, Fernandez-Gamiz U, Noeiaghdam S, Shah SA, et al. Numerical analysis of unsteady hybrid nanofluid flow comprising CNTs-ferrousoxide/water with variable magnetic field. Nanomaterials. 2022;12(2):180. 10.3390/nano12020180Search in Google Scholar PubMed PubMed Central

[15] Seydel R, Hlavacek V. Role of continuation in engineering analysis. Chem Eng Sci. 1987;42(6):1281–95. 10.1016/0009-2509(87)85001-7Search in Google Scholar

[16] Mkhatshwa M, Khumalo M. Irreversibility scrutinization on EMHD Darcy-Forchheimer slip flow of Carreau hybrid nanofluid through a stretchable surface in porous medium with temperature-variant properties. Heat Transfer. 2023;52(1):395–429. 10.1002/htj.22700Search in Google Scholar

[17] Mkhatshwa MP. Overlapping grid spectral collocation methods for nonlinear differential equations modelling fluid flow problems. [Ph.D. thesis 2020]. University of KwaZulu-Natal; 2023. Search in Google Scholar

[18] Mkhatshwa M, Motsa S, Sibanda P. Overlapping multi-domain bivariate spectral method for systems of nonlinear PDEs with fluid mechanics applications. In: Advances in Fluid Dynamics: Selected Proceedings of ICAFD 2018. Springer; 2021. p. 685–99. 10.1007/978-981-15-4308-1_54Search in Google Scholar

[19] Mkhatshwa M. Overlapping grid spectral collocation approach for electrical MHD bioconvection Darcy-Forchheimer flow of a Carreau-Yasuda nanoliquid over a periodically accelerating surface. Heat Transfer. 2022;51(2):1468–500. 10.1002/htj.22360Search in Google Scholar

[20] Mkhatshwa MP, Motsa SS, Sibanda P. Numerical solution of time-dependent Emden-Fowler Eqs. using bivariate spectral collocation method on overlapping grids. Nonlinear Eng. 2020;9(1):299–318. 10.1515/nleng-2020-0017Search in Google Scholar

[21] Motsa S. On the optimal auxiliary linear operator for the spectral homotopy analysis method solution of nonlinear ordinary differential equations. Math Problems Eng. 2014;2014:697845. 10.1155/2014/697845Search in Google Scholar

[22] Motsa S, Magagula V, Makukula Z. Simple iteration methods for non-linear differential equations: Theory and Development. 10th Annual Research Workshop on Numerical Methods for Differential Equations, University of KwaZulu-Natal, KwaZulu-Natal, South Africa. 2017. Search in Google Scholar

[23] Motsa SS. A new spectral local linearization method for nonlinear boundary layer flow problems. J Appl Math. 2013;2013. 10.1155/2013/423628. Search in Google Scholar

[24] Otegbeye O, Ansari MS. A finite difference based simple iteration method for solving boundary layer flow problems. In: AIP Conference Proceedings. vol. 2435. AIP Publishing; 2022. 10.1063/5.0084396Search in Google Scholar

[25] Hayat T, Qasim M, Abbas Z. Homotopy solution for the unsteady three-dimensional MHD flow and mass transfer in a porous space. Commun Nonlinear Sci Numer Simulat. 2010;15(9):2375–87. 10.1016/j.cnsns.2009.09.013Search in Google Scholar

[26] Abbas Z, Javed T, Sajid M, Ali N. Unsteady MHD flow and heat transfer on a stretching sheet in a rotating fluid. J Taiwan Instit Chem Eng. 2010;41(6):644–50. 10.1016/j.jtice.2010.02.002Search in Google Scholar

[27] Ahmed SE, Hussein AK, Mansour M, Raizah ZA, Zhang X. MHD mixed convection in trapezoidal enclosures filled with micropolar nanofluids. Nanosci Tech Int J. 2018;9(4):343–72. 10.1615/NanoSciTechnolIntJ.2018026118Search in Google Scholar

[28] Lund LA, Omar Z, Khan I, Sherif ESM. Dual branches of MHD three-dimensional rotating flow of hybrid nanofluid on nonlinear shrinking sheet. Comput Materials Contin. 2020;66(1):127–39. 10.32604/cmc.2020.013120Search in Google Scholar

[29] Suresh S, Venkitaraj K, Selvakumar P, Chandrasekar M. Synthesis of Al2O3–Cu/water hybrid nanofluids using two step method and its thermo physical properties. Colloids Surfaces A Physicochem Eng Aspects. 2011;388(1–3):41–8. 10.1016/j.colsurfa.2011.08.005Search in Google Scholar

[30] Madhukesh JK, Ramesh GK, Roopa GS, Prasannakumara BC, Shah NA, Yook SJ. 3D flow of hybrid nanomaterial through a circular cylinder: Saddle and Nodal point aspects. Mathematics. 2022;10(7):1185. 10.3390/math10071185Search in Google Scholar

[31] Rana S, Nawaz M, Alaoui M. Three dimensional heat transfer in the Carreau-Yasuda hybrid nanofluid with Hall and ion slip effects. Phys Scr. 2021;96(12):125215. 10.1088/1402-4896/ac2379Search in Google Scholar

[32] Williams III J, Rhyne T. Boundary layer development on a wedge impulsively set into motion. SIAM J Appl Math. 1980;38(2):215–24. 10.1137/0138019Search in Google Scholar

[33] Magagula VM, Motsa SS, Sibanda P, Dlamini PG. On a bivariate spectral relaxation method for unsteady magneto-hydrodynamic flow in porous media. SpringerPlus. 2016;5(1):1–15. 10.1186/s40064-016-2053-4Search in Google Scholar PubMed PubMed Central

[34] Revers M. On the approximation of certain functions by interpolating polynomials. Bulletin Aust Math Soc. 1998;58(3):505–12. 10.1017/S0004972700032494Search in Google Scholar

[35] Devi SU, Devi SA. Heat transfer enhancement of Cu-Al2O3/water hybrid nanofluid flow over a stretching sheet. J Nigerian Math Soc. 2017;36(2):419–33. Search in Google Scholar

[36] Khashi’ie NS, Arifin NM, Hafidzuddin EH, Wahi N. Thermally stratified flow of Cu-Al2O3/water hybrid nanofluid past a permeable stretching/shrinking circular cylinder. J Adv Res Fluid Mech Therm Sci. 2019;63(1):154–63. Search in Google Scholar

[37] Waini I, Ishak A, Pop I. Hybrid nanofluid flow and heat transfer past a permeable stretching/shrinking surface with a convective boundary condition. In: Journal of Physics: Conference Series. vol. 1366. IOP Publishing; 2019. p. 012022. 10.1088/1742-6596/1366/1/012022Search in Google Scholar

Received: 2023-06-12

Revised: 2023-11-28

Accepted: 2023-12-01

Published Online: 2025-09-02

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

Articles in the same Issue

https://doi.org/10.1515/nleng-2022-0368

Keywords for this article

hybrid nanofluid; multi-domain overlapping grid; bivariate spectral simple iteration method

Creative Commons

BY 4.0