Nowadays, researchers have created novel fluids to suit the requirement for improved heat transmission and thermal conductivity. The heat transfer properties of conventional fluids like water, ethylene glycol, glycerin, and ethanol are limited in applications such as power generation, chemical processes, and heating and cooling systems. Scientists are investigating better heat transmission materials to meet the rising energy demands and address concerns about shortages of resources and environmental impacts. A single kind of nanoparticle can be added to the aforementioned fluids to make up for their deficiency. The term nanofluid refers to this process, which was initially studied by Choi and Eastman. ¹ Researchers have considered many combinations of nanoparticles, including semiconductors ( SiO ₂ , TiO ₂), metallic oxides ( Al ₂ O ₃ , CuO), and metal nanoparticles ( Al, Cu, and Fe). Heat transfer through nanofluids has been investigated from various perspectives. Mahian et al. ² presented the uses of nanofluid in a variety of contexts, such as renewable energy systems. They also covered the advantages of energy systems from an environmental perspective when employing nanofluid. Mansoury et al. ³ investigated the flow of Al ₂ O ₃ /H ₂ O nanofluids through parallel heat exchangers. Unfortunately, to obtain the required thermal performance, a single nanoparticle suspension is insufficient. To attain the required thermal properties, hybrid nanofluid is employed. Several experimental and theoretical models have been published and studied in order to use hybrid nanofluids for more efficient industrial and technological processes. Mahanthesh et al. ⁴ reviewed the flow behavior of hybrid nanofluids, focusing on the effects of Brownian motion and thermophoresis. Sensitivity analysis revealed that the Brownian motion parameter has the most significant impact on the heat transfer rate. Bilal et al. ⁵ investigated the Darcy-Forchheimer mixed convection flow of hybrid nanofluids through an inclined, extending cylinder using the homotopy analysis method. Their findings indicated that CNT−Fe ₃O ₄/H ₂O hybrid nanofluids enhance the thermal efficiency of the base fluid more effectively than conventional fluids. The investigation provides further information on the hybrid nanofluid flow as provided in these reference. ^{6, 7}

Scientists and engineers have been motivated by the non-Newtonian motion of fluids because these materials have multiple applications in science and technology processes. Some examples of non-Newtonian fluid are Paints, mud, Soap, glues, apple sauce, printing ink, shampoos, sugar solutions, tomato paste, etc. Furthermore, non-Newtonian fluids are used in several kinds of fields of study, such as the chemical engineering field, biology, and geophysical sciences. The Williamson model is the fluid model that is being studied. In the case of Williamson fluid rheology, the constitutive equation for stress-strain is non-linear. Therefore, viscosity that depends on the shear rate forms the basis of the Williamson fluid model. This model captures both shear-thinning and shear-thickening behavior. Williamson ⁸ introduced the Williamson fluid model, which defines the relationship between stress and strain in pseudoplastic materials. Kebede et al. ⁹ studied the heat and mass transport properties of Williamson nanofluid flow. Kumar et al. ¹⁰ explored the impact of Brownian motion and thermophoresis on heat and mass transfer in pseudoplastic materials. Reddy et al. ¹¹ investigated the impact of Cattaneo–Christov heat flux and suction/injection on the entropy generation in the system. Chakradhar et al. ¹² examined the influence of MHD effects on the peristaltic motion of Williamson fluid within a porous channel.

Heat and mass transport are natural phenomena caused by concentration and temperature differences within or between materials. Many industrial processes, including wire drawing, artificial fiber production, paper manufacturing, chemical waste migration, and distillation, are influenced by this phenomenon. Significant efforts have been made in the past to study heat and mass transport mechanisms using Fourier’s law of heat transfer ¹³ and Fick’s law of diffusion. ¹⁴ Previously,”Fourier’s law of heat conduction” was commonly used to explain heat transfer. However, this law does not fully capture the fundamental nature of heat transfer. To address this issue, Cattaneo ¹⁵ introduced thermal relaxation into Fourier’s theory, making heat transfer resemble thermal wave propagation at normal speed. Christov ¹⁶ improved the Cattaneo model’s thermal relaxation time by using upper-convected Oldroyd derivatives for frame-invariant formation. Sui et al. ¹⁷ extended the Cattaneo-Christov model to mass diffusion problems, applying it to the Maxwell nanofluid mass diffusion across a moving surface. Vinodkumar et al. ¹⁸ explored the behavior of MHD convective flow of non-Newtonian nanofluids over a chemically reacting porous sheet. Vinodkumar et al. ¹⁹ examined the effects of higher-order chemical reactions and radiation on the MHD flow of a Maxwell nanofluid. Recent developments on this concept are gathered in. ^{20– 25}

Magnetohydrodynamics is the study of the interaction between a magnetic field and electrically conducting fluids, such as salt water solutions, liquid metals, and plasmas. Magnetohydrodynamic flow studies have significant applications in chemistry, physics, and engineering. Additionally, MHD flow has applications in industrial equipment, pumps, electric transformers, and hydro-magnetic generators, among others. Zhao et al. ²⁶ studied the thermally induced electro-kinetic flow of Al ₂O ₃-water nanofluid through a permeable micro-tube to examine its heat transfer properties, considering the effect of an applied magnetic field. Chamkaa et al. ²⁷ used the finite-difference method to analyze mixed convection in a square cavity filled with Cu-water nanofluid. MS et al. ²⁸ analyzed the entropy optimization in the MHD flow of non-Newtonian nanofluids with chemical reactions and thermal energy effects. Vinodkumar et al. ²⁹ investigated the activation and behavior of energy in the magnetohydrodynamic (MHD) flow of Casson nanofluids through a porous medium. Recent studies on MHD can be found in references. ^{30, 31}

Flow over a stretching sheet is a classic fluid mechanics problem with various practical applications in industries such as glass fiber production, glass blowing, wire drawing, and copper wire tinning. The concept of boundary layer flow for solid surface motion at a constant speed was first introduced by Sakiadis. ^{32– 34} Several researchers have recently investigated the flow of boundary layers over stretchable surfaces. ^{35, 36}

In fluid dynamics, a porous medium refers to a solid material that contains a network of interconnected voids or pores, allowing fluids (such as gases or liquids) to flow through it. The structure of the porous medium, including the size, shape, and distribution of its pores, significantly influences the flow behavior and other transport phenomena (e.g., heat and mass transfer) within the medium. Shaw et al. ³⁷ investigated the effects of entropy production in a Casson fluid containing MWCNT/Fe ₃ O ₄ through a stretched disc and the Darcy-Forchheimer porous medium concept using a numerical approach. There are further investigations on the porous media in. ^{38, 39}

Nonlinear thermal radiation refers to the heat transfer process where the radiation heat flux does not follow a simple linear relationship with temperature. Unlike linear radiation, where the heat flux is proportional to the fourth power of temperature (as in the Stefan-Boltzmann law), nonlinear thermal radiation accounts for more complex interactions between the radiating bodies and the surrounding medium. This phenomenon arises in high-temperature environments, where the radiation intensity varies non-linearly due to factors such as temperature gradients, optical properties of materials, or the nature of the medium. Alamirew et al. ⁴⁰ studied the impact of nonlinear thermal radiation, ion slip, and Hall on MHD Williamson nanofluid flow over a stretched sheet. Building on this, other researchers ^{41, 42} have explored the effects of thermal radiation in various geometries and methods.

Several studies have investigated Williamson nanofluid flow over a stretching sheet using various geometries and methods. Different approaches have been employed in the literature to study the heat and mass fluxes in Williamson nanofluid flow. However, to the best of the authors’ knowledge, no studies have yet considered the combined effects of non-linear thermal radiation, Joule heating, non-Fourier heat flux, non-Fick mass flux, chemical reactions, heat generation/absorption, Brownian motion, and thermophoresis. The present study aims to provide a comprehensive analysis of mass and heat transfer in 3D MHD Williamson hybrid nanofluid flow, which consists of Cu and Al ₂O ₃ nanoparticles in ethylene glycol, over a linearly stretching porous sheet. By applying similarity transformations, the governing nonlinear partial differential equations are converted into a system of ordinary differential equations. The reduced mathematical model is then solved using the bvp4c function in MATLAB. The study uses tables and figures to analyze the effects of various parameters, such as velocity, temperature, concentration, skin friction, the Nusselt number, and the Sherwood number. The findings of this study have important implications for engineering and industrial applications, including aerospace engineering, nuclear power plants, electronic cooling, biomedical systems, and the polymer and textile industries.

Assume a steady three-dimensional hybrid nanofluid ( Cu − Al 2 O 3 /Ethylene glycol) flow past a stretching porous sheet in the presence of an applied magnetic field. The schematic representation in Cartesian coordinates is displayed in Figure 1. The flow is incompressible and laminar. Cattaneo-Christov Heat and Mass Flux Model, Non-Linear Thermal Radiation, Joule Heating, Viscous Dissipation, Porous Medium and Magnetic Field Effects, heat generation, and chemical reaction are applied. The velocities of the stretching sheets are v w = by in y and u w = ax in x , and a constant external magnetic field B is applied in z . The temperature and concentration of the surface are kept constant at T w and C w respectively, which is higher than its ambient temperature, T ∞ and ambient concentration, C ∞ . Furthermore, the double diffusion Cattaneo-Christove theorya generalization of the well-known Fourier and Fick laws was used to study heat and mass transmission. ^{43, 44} q → + κ ∇ T + λ e [ q → t + V → · ∇ q → − q → · ∇ V → + ( ∇ · V → ) q → ] = 0 (1) J → + D B ∇ C + λ c [ J → t + V → · ∇ J → − J → · ∇ V → + ( ∇ · V → ) J → ] = 0 (2)

In equation (1) and (2), q → and J → , D B , T , C , κ , V → , λ e and λ c represents the heat and mass fluxes, the Brownian motion, the temperature, the concentration, the thermal conductivity, the velocity field, the energy and concentration relaxation factors, respectively. By letting λ e = 0 = λ c , the famous Fourier’s and Fick’s laws can be obtained. From equations (1) and (2), which generalize the classic Fourier and Fick’s laws, the steady-state conditions of the nanofluid flow are expressed as follows: q → + κ ∇ T + λ e [ V → · ∇ q → − q → · ∇ V → + ( ∇ · V → ) q → ] = 0 (3) J → + D B ∇ C + λ c [ V → · ∇ J → − J → · ∇ V → + ( ∇ · V → ) J → ] = 0 . (4)

The Williamson fluid model is given Refs. 45, 46 by τ = − P → I + S → (5) S → = [ μ ∞ + μ 0 − μ ∞ 1 − Γ γ ̇ ] A 1 (6)

The shear rate ( γ ̇ ) and is defined as γ ̇ = 0.5 trace ( A 1 ) 2 (7)

This is taken into consideration for pseudo-plastic fluids μ ∞ = 0 and Γ γ ̇ < 1 . Equation (6) can be written as S → = [ μ 0 / ( 1 − Γ γ ̇ ) ] A 1 (8)

Using the binomial expansion to Equation (8), given as S → ≈ μ 0 [ 1 + Γ γ ̇ ] A 1 . (9)

Considering the above assumptions, and after applying the boundary layer approximations, the governing equations are expressed as follows ^{43, 47, 48} ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z = 0 (10) u ∂ u ∂ x + v ∂ u ∂ y + w ∂ u ∂ z = μ hnf ρ hnf ∂ 2 u ∂ z 2 + 2 Γ μ hnf ρ hnf ∂ u ∂ z ∂ 2 u ∂ z 2 − σ hnf B 2 ρ hnf u − μ hnf ρ hnf K p u (11) u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z = μ hnf ρ hnf ∂ 2 v ∂ z 2 + 2 Γ μ hnf ρ hnf ∂ v ∂ z ∂ 2 v ∂ z 2 − σ hnf B 2 ρ hnf v − μ hnf ρ hnf K p v (12) u ∂ T ∂ x + v ∂ T ∂ y + w ∂ T ∂ z = α hnf ∂ 2 T ∂ z 2 − 1 ( ρCp ) hnf ∂ q r ∂ z − λ e [ u 2 ∂ 2 T ∂ x 2 + v 2 ∂ 2 T ∂ y 2 + w 2 ∂ 2 T ∂ z 2 + 2 ( uv ∂ 2 T ∂ x ∂ y + vw ∂ 2 T ∂ y ∂ z + uw ∂ 2 T ∂ x ∂ z ) + ( u ∂ u ∂ x + v ∂ u ∂ y + w ∂ u ∂ z ) ∂ T ∂ x + ( u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z ) ∂ T ∂ y + ( u ∂ w ∂ x + v ∂ w ∂ y + w ∂ w ∂ z ) ∂ T ∂ z ] + Q 0 ( T − T ∞ ) ( ρCp ) hnf + σ hnf B 2 ( u 2 + v 2 ) ( ρCp ) hnf + τ [ D B ∂ C ∂ z ∂ T ∂ z + D T T ∞ ( ∂ T ∂ z ) 2 ] (13) u ∂ C ∂ x + v ∂ C ∂ y + w ∂ C ∂ z + κ c λ c = D B ∂ 2 C ∂ z 2 + D T T ∞ ( ∂ 2 T ∂ z 2 ) − K r ( C − C ∞ ) − λ c [ u 2 ∂ 2 C ∂ x 2 + v 2 ∂ 2 C ∂ y 2 + w 2 ∂ 2 C ∂ z 2 + 2 ( uv ∂ 2 C ∂ x ∂ y + vw ∂ 2 C ∂ y ∂ z + uw ∂ 2 C ∂ x ∂ z ) + ( u ∂ u ∂ x + v ∂ u ∂ y + w ∂ u ∂ z ) ∂ C ∂ x + ( u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z ) ∂ C ∂ y + ( u ∂ w ∂ x + v ∂ w ∂ y + w ∂ w ∂ z ) ∂ C ∂ z ] (14)

The boundary conditions for the problem are given by ^{47, 49} { w = 0 , v = v w = by , u = u w = ax , T = T w , C = C w , at z → 0 u = v = 0 , T = T ∞ , C = C ∞ , as z → ∞ (15)

The hybrid nanofluids thermophysical characteristics are described as ^{46, 47} : μ hnf = μ f ( 1 − ϕ 1 ) 2.5 ( 1 − ϕ 2 ) 2.5 , α hnf = κ hnf ( ρCp ) hnf , ρ hnf = ( 1 − ϕ 2 ) [ ( 1 − ϕ 1 ) ρ f + ϕ 1 ρ n 1 ] + ϕ 2 ρ n 2 , ( ρCp ) hnf = ( 1 − ϕ 2 ) [ ( 1 − ϕ 1 ) ( ρCp ) f + ϕ 1 ( ρCp ) n 1 ] + ϕ 2 ( ρCp ) n 2 , κ hnf κ nf = κ n 2 + 2 κ nf − 2 ϕ 2 ( κ nf − κ n 2 ) κ n 2 + 2 κ nf + ϕ 2 ( κ nf − κ n 2 ) , κ nf κ f = κ n 1 + 2 κ f − 2 ϕ 1 ( κ f − κ n 1 ) κ n + 2 κ f + ϕ ( κ f − κ n 1 ) , σ hnf σ nf = σ n 2 + 2 σ nf − 2 ϕ 2 ( σ nf − σ n 2 ) σ n 2 + 2 σ nf + ϕ 2 ( σ nf − σ n 2 ) , σ nf σ f = σ n 1 + 2 σ f − 2 ϕ 1 ( σ f − σ n 1 ) σ n + 2 σ f + ϕ ( σ f − σ n 1 ) , where μ f , ϕ 1 , ϕ 2 , ρ f , κ , ( ρCp ) f , σ hnf , hnf , nf , f , and n the base fluid viscosity, nanoparticle volume fraction of Cu, nanoparticle volume fraction of Al 2 O 3 , density, thermal conductivity, heat capacitance, electrical conductivity of hybrid nanofluid, hybrid nanofluid, nanofluid, base fluid, and nanoparticle. Furthermore, in terms of shape, only spherical nanoparticles are taken into consideration. This hybrid nanofluid is applicable in aerospace, nuclear reactors, biomedical engineering, electronic cooling, and renewable energy systems, making it an optimal choice for advanced heat transfer applications.

Using the nonlinear Rosseland diffusion approximation, the radiative heat flux q r can be expressed as: ⁴⁰ q r = − 4 σ ∗ 3 κ ∗ ∂ T 4 ∂ z (16) q r = − 16 σ ∗ 3 κ ∗ T 3 ∂ T ∂ z (17) and therefore ∂ q r ∂ z = − 16 σ ∗ 3 κ ∗ ∂ ∂ z ( T 3 ∂ T ∂ z ) (18)

By means of Equations (16), (17) and (18), Equation (13) reduces to u ∂ T ∂ x + v ∂ T ∂ y + w ∂ T ∂ z = α hnf ∂ 2 T ∂ z 2 + 1 ( ρCp ) hnf 16 σ ∗ 3 κ ∗ ∂ ∂ z ( T 3 ∂ T ∂ z ) − λ e [ u 2 ∂ 2 T ∂ x 2 + v 2 ∂ 2 T ∂ y 2 + w 2 ∂ 2 T ∂ z 2 + 2 ( uv ∂ 2 T ∂ x ∂ y + vw ∂ 2 T ∂ y ∂ z + uw ∂ 2 T ∂ x ∂ z ) + ( u ∂ u ∂ x + v ∂ u ∂ y + w ∂ u ∂ z ) ∂ T ∂ x + ( u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z ) ∂ T ∂ y + ( u ∂ w ∂ x + v ∂ w ∂ y + w ∂ w ∂ z ) ∂ T ∂ z ] + Q 0 ( T − T ∞ ) ( ρCp ) hnf + σ hnf B 2 ( u 2 + v 2 ) ( ρCp ) hnf + τ [ D B ∂ C ∂ z ∂ T ∂ z + D T T ∞ ( ∂ 2 T ∂ z 2 ) 2 ] (19)

The similarity transformations are given by ^{47, 49} { w = − a ν f ( f ( η ) + g ( η ) ) , u = ax f ′ ( η ) , v = ay g ′ ( η ) , η = z a ν f T = T ∞ + ( T w − T ∞ ) θ ( η ) , C = C ∞ + ( C w − C ∞ ) Φ ( η ) (20)

Equations (10) to (14) are simplified accordingly: A 1 A 2 f ′ ′ ′ [ 1 + We f ′ ′ ] + ( f + g ) f ′ ′ − f ′ 2 − A 3 A 2 M f ′ − A 1 A 2 K f ′ = 0 (21) A 1 A 2 g ′ ′ ′ [ 1 + We g ′ ′ ] + ( f + g ) g ′ ′ − g ′ 2 − A 3 A 2 M g ′ − A 1 A 2 K g ′ = 0 (22) A 4 A 5 Pr θ ′ ′ + ( f + g ) θ ′ − β 1 [ ( f + g ) 2 θ ′ ′ + ( f + g ) ( f ′ + g ′ ) θ ′ ] + 1 A 5 Pr R [ ( 1 + ( θ w − 1 ) ) 2 × 3 ( θ w − 1 ) θ ′ 2 + ( 1 + ( θ w − 1 ) θ ) 3 θ ′ ′ ] + 1 A 5 Qθ + Nb θ ′ Φ ′ + Nt θ 2 + A 3 A 5 M ( Ecx ( f ′ ) 2 + Ecy ( g ′ ) 2 ) = 0 (23) Φ ′ ′ + Nt Nb θ ′ ′ + Sc ( f + g ) Φ ′ − β 2 Sc [ ( f + g ) 2 Φ ′ ′ + ( f + g ) ( f ′ + g ′ ) Φ ′ ] − ScH Φ = 0 (24) with subject to boundary conditions { f ′ ( 0 ) = 1 , f ( 0 ) = 0 , g ( 0 ) = 0 , g ′ ( 0 ) = d , θ ( 0 ) = 1 , Φ ( 0 ) = 1 , at η = 0 f ′ ( η ) = 0 , g ′ ( η ) = 0 , θ ( η ) = 0 , Φ ( η ) = 0 , as η → ∞ (25) where, the Weissenberg number ( We ) , magnetic field parameter ( M ) , porosity parameter ( K ) , Prandtl number ( Pr ) , non-linear thermal radiation ( R ) , ratio of temperature ( θ w ) , Q is the heat generation ( Q > 0 ) or absorption parameter ( Q < 0 ) , Brownian motion parameter ( Nb ) , Thermophoresis parameter ( Nt ) , Eckert number ( Ec ) along the x - and y -axis, thermal relaxation parameter ( β 1 ) , concentration relaxation parameter ( β 2 ) , Schmidt number ( Sc ) , Chemical reaction ( H ) , d is ratio of stretching sheet, and A 1 , A 2 , A 3 , A 4 , and A 5 are dynamic viscosity, density, electrical conductivity, thermal conductivity, and heat capacitance of the hybrid nanofluid, respectively can be described as follows: We = x Γ 2 a 3 ν f , M = σ f B 2 ρ f a , K = ν f a K p , Pr = ν f ( ρCp ) f κ f , R = 16 σ ∗ T ∞ 3 3 κ κ ∗ , θ w = T w T ∞ Q = Q 0 a ( ρCp ) f , Nb = D B τ ( C w − C ∞ ) ν f , Nt = D T τ ( T w − T ∞ ) ν f T ∞ , E c x = u w 2 ( T w − T ∞ ) ( Cp ) f E c y = v w 2 ( T w − T ∞ ) ( C p ) f , β 1 = λ e a , β 2 = λ c a , Sc = ν f D B , H = K r a , d = b a A 1 = μ hnf μ f , A 2 = ρ hnf ρ f , A 3 = σ hnf σ f , A 4 = κ hnf κ f , A 5 = ( ρCp ) hnf ( ρCp ) f

The engineering components of the skin friction coefficients C fx and C fy , Nusselt number N u x and Sherwood number S h x are defined as C fx = τ xz ρ f u w 2 , C fy = τ yz ρ f v w 2 , N u x = x q w κ f ( T w − T ∞ ) , S h x = x j m D B ( C w − C ∞ ) , (26) where τ xz and τ yz are the wall shear stresses in the x- and y-direction respectively, given as τ xz = μ hnf [ ∂ u ∂ z + Γ 2 ( ∂ u ∂ z ) 2 ] z = 0 , τ yz = μ hnf [ ∂ v ∂ z + Γ 2 ( ∂ v ∂ z ) 2 ] z = 0 (27) and the wall heat and mass flux from the sheets, which is given by: q w = − [ κ hnf + 16 σ ∗ T ∞ 3 3 κ ∗ ] ( ∂ T ∂ z ) z = 0 , j m = − D B ( ∂ C ∂ z ) z = 0 (28)

By using Equations (26), (27) and (28), the dimensionless variables, we obtain: ( R e x ) 1 2 C fx = A 1 [ 1 + We 2 f ′ ′ ( 0 ) ] f ′ ′ ( 0 ) (29) d 1.5 ( R e y ) 1 2 C fy = A 1 [ 1 + We 2 g ′ ′ ( 0 ) ] g ′ ′ ( 0 ) (30) ( R e y ) − 1 2 N u x = − ( A 4 + R ) θ ′ ( 0 ) (31) ( R e y ) − 1 2 S h x = − Φ ′ ( 0 ) (32) where R e y = y v w ν f and R e x = x u w ν f are the local Reynolds numbers.

Using Table 1, and the MATLAB software ( https://github.com/asfawmat/BVP-MATLAB-Implementation) uses the bvp4c technique to numerically solve equations (21) through (24) and the boundary condition (25). The bvp4c is an effective instrument that provides precision and dependability in addressing boundary value problems for differential equation systems. BVP4C is chosen for its high accuracy, adaptive mesh refinement, and robustness in handling boundary value problems (BVPs), making it well-suited for nonlinear and coupled ODEs arising in Williamson hybrid nanofluid flow. Unlike shooting methods, BVP4C provides stability and effectively handles multi-point boundary conditions, making it ideal for fluid flow and heat transfer problems. Although methods like Runge-Kutta and finite differences were considered, they proved less efficient due to stability and convergence issues. However, BVP4C requires careful selection of initial guesses due to its computational intensity. With a high degree of accuracy specifically, fourth- order accuracy this solver exemplifies a collocation technique that provides a precise, continuous solution. The mesh selection and error control are determined by evaluating the residual of the continuous solution. As discussed later, the boundary value problems must be transformed into a system of first-order initial value problems (IVPs) to apply the bvp4c method using the shooting technique. Now let us defined the new variable by the equation { y 1 = f , y 2 = f ′ , y 3 = f ′ ′ , y 4 = g , y 5 = g ′ , y 6 = g ′ ′ , y 7 = θ , y 8 = θ ′ , y 9 = Φ , y 10 = Φ ′ . (33) y 1 ′ = y 2 y 2 ′ = y 3 y 3 ′ = 1 A 1 A 2 ( 1 + We y 3 ) [ A 3 A 2 M y 2 + A 1 A 2 K y 2 + y 2 2 − ( y 1 + y 4 ) y 3 ] , y 4 ′ = y 5 , y 5 ′ = y 6 , y 6 ′ = 1 A 1 A 2 ( 1 + We y 6 ) [ A 3 A 2 M y 5 + A 1 A 2 K y 5 + y 5 2 − ( y 1 + y 4 ) y 6 ] , y 7 ′ = y 8 , y 8 ′ = − [ 1 ( A 4 A 5 Pr + R A 5 Pr ( 1 + ( θ w − 1 ) y 7 ) 3 − β 1 ( y 1 + y 4 ) 2 ) ] × [ R A 5 Pr ( 1 + ( θ w − 1 ) ) 2 × 3 ( θ w − 1 ) y 8 2 + ( y 1 + y 4 ) y 8 − β 1 ( y 1 + y 4 ) ( y 2 + y 5 ) y 8 + 1 A 5 Q y 7 + Nb y 8 y 10 + Nt y 8 2 + A 3 A 5 M ( E c x ( y 2 ) 2 + E c y ( y 5 ) 2 ) ] , y ′ 9 = y 10 , y 10 ′ = [ 1 ( 1 − β 2 Sc ( y 1 + y 4 ) 2 ) ] × [ β 2 Sc ( ( y 1 + y 4 ) ( y 2 + y 5 ) − Sc ( y 1 + y 4 ) ) y 10 + ScH y 9 − Nt Nb [ − [ 1 ( A 4 A 5 Pr + R A 5 Pr ( 1 + ( θ w − 1 ) θ ) 3 − β 1 ( y 1 + y 4 ) 2 ) ] × [ R A 5 Pr ( ( 1 + ( θ w − 1 ) ) 2 × 3 ( θ w − 1 ) y 8 2 ) + ( y 1 + y 4 ) y 8 − β 1 ( y 1 + y 4 ) ( y 2 + y 5 ) y 8 + 1 A 5 Q y 7 + Nb y 8 y 10 + Nt y 8 2 + A 3 A 5 M ( E c x ( y 2 ) 2 + E c y ( y 5 ) 2 ) ] ] with corresponding initial conditions { y 1 ( 0 ) = 0 , y 2 ( 0 ) = 1 , y 3 ( 0 ) = α 1 , y 4 ( 0 ) = 0 , y 5 ( 0 ) = d y 6 ( 0 ) = α 2 , y 7 ( 0 ) = 1 , y 8 ( 0 ) = α 3 , y 9 ( 0 ) = 1 , y 10 ( 0 ) = α 4 (34) where α 1 , α 2 , α 3 , and α 4 are the missing initial conditions.

Utilizing the bvp4c technique, a system of differential equations of the form y ′ = f ( x , y ) is integrated according to the specified boundary conditions. The bvp4c routine uses finite difference method with an achievable accuracy of about 10 − 7 .

In this section, we present the numerical solutions to the problem, considering various physical effects. The system of equations (21)- (24) is solved using the numerical bvp4c method, which satisfies the boundary conditions given by equation (25). Additionally, we conducted a comparative analysis of our numerical results with those from previous studies to verify the accuracy of our solution and evaluate its consistency with earlier findings. The comparison of our findings with those from earlier research is presented in Table 2. Comparing our findings with those of previous studies, the comparison’s results show generally excellent agreement. Using a variety of graphs, the impacts of various physical factors on temperature, concentration, velocity, mass transfer rate, surface drag coefficient, and the hybrid Cu-Al ₂ O ₃ /ethylene glycol nanofluid phase and nanofluid Al ₂ O ₃ /ethylene glycol phase are displayed. In each figure, a comparison of mono- and hybrid nanofluids is also provided. The numerical computations are discussed by keeping d = 0.5 , Sc = 1 , Nb = Nt = 0.1 , We = 0.2 , R = θ w = 1.2 , K = 0.1 , Pr = 1 , Ec = 0.2 , β 1 = β 2 = 0.1 , H = 0.2 , Q = 0.1 , M = 1,0.01 ≤ ϕ 1 ≤ 0.05 , 0.01 ≤ ϕ 2 ≤ 0.05 throughout the complete study.

An analysis of the differences between the numerical results of Wang, ⁵¹ You & Wang, ⁴⁷ and the current study’s conclusions by calculating the values of f " ( 0 ) and g " ( 0 ) under the conditions of Pr = 1, Nb = 1 × 10 −14, when We = M = K = β 1 = β 2 = R = θ w = Q = Nt = Ec = Sc = H = ϕ 1 = ϕ 2 = 0, as shown in Table 2 below.

This section presents the results for a 3D MHD Williamson hybrid nanofluid flow. The study utilizes the Cattaneo-Christov and modified Boungerno’s models to simulate the Cu−Al ₂O ₃/Ethylene glycol MHD Williamson hybrid nanofluid flow over a linearly stretching sheet, influenced by factors such as magnetic field, heat generation/absorption, non-linear thermal radiation, Joule heating, Brownian motion, thermophoresis, porosity, Williamson fluid parameter, Eckert number, Prandtl number, Schmidt number, chemical reaction parameter, stretching ratio, thermal relaxation, and concentration relaxation parameters. By applying appropriate similarity variables, the system of partial differential equations is transformed into a non-linear ordinary differential equation formulation. These ordinary differential equations are solved numerically using the bvp4c solver on the MATLAB platform. Similarly, graphical analysis is used to plot the obtained flow parameters so that their effects on flow, heat, and mass can be easily seen. Researchers and engineers working on related issues can benefit greatly from the obtained data. Moreover, our findings demonstrate the effectiveness and utility of hybrid nanofluid performance.

The present study leads to the following deductions: •

When the stretching ratio parameter rises, the primary velocity profile diminishes, but the secondary velocity profile grows.

A dedicated section on future research directions follows the conclusion, outlining potential extensions of this study. These include examining unsteady flow conditions to capture transient effects, exploring various hybrid nanoparticle combinations to enhance heat transfer performance, and extending the model to incorporate non-Newtonian rheological properties relevant to biomedical and industrial applications.

Asfaw Tsegaye Moltot contributed to Writing review & editing, conceptualization, methodology, formal analysis, validation, and writing the original draft. Eshetu Haile contributed to Writing the review, editing, supervision, conceptualization, and resources. Gurju Awgichew contributed to Writing the review, editing, supervision and resources. Hunegnaw Dessie contributed to Writing the review, editing, supervision and resources.

Ethical approval and consent were not required.