1 Introduction

Transport procedures through a porous space with high Reynolds numbers are commonly encountered in various chemical, mechanical, geophysical, electrochemical, and metallurgical processes. Initially, Darcy^[1] proposed a law to explicate the dynamic phenomenon in porous media by examining the flow of sandy water through pebbles by considering a low Reynolds number. Later, he worked to develop a similar flow description for high Reynolds numbers, but was unsuccessful in this attempt. After studying and analyzing the pioneering work done by Darcy, Forchheimer^[2] overcame the shortcoming in the Darcy law by including the inertial term in the momentum equation. This remarkable improvement in the law opened a new area of research in porous media and set a benchmark for further studies in this field. The additive term in the Darcy law was renowned by Muskat^[3] as the Forchheimer term. The Darcy-Forchheimer model within the frame of a stretched linearized configuration in the presence of an electric field was considered by Pal and Mondal^[4]. They concluded that the electric field tends to decrease the nanoparticle volume fraction. Ganesh et al.^[5] executed the aspects of magneto nanomaterial in an absorbent surface. They implemented a modified Forchheimer relation to investigate the porosity impact on the fluid flow. Haq et al.^[6] explored the natural convection phenomenon for magnetohydrodynamic (MHD) flow in a porous corrugated cavity at various wavelengths and a partially heated domain by the finite element method (FEM). They presented the flow patterns and the temperature distribution within the entire domain of the cavity through streamlines and isotherms. Hayat et al.^[7] elucidated an improved Forchheimer impression on the Maxwell fluid exposed to a Cattaneo-Christov expression. In recent years, the Darcy-Forchheimer flow in various physical conditions has been interpreted by researchers. Muhammad et al.^[8] addressed the flow behavior of Maxwell nanofluid in porous media by implementing a zero-mass flux condition. They leveraged the Darcy-Forchheimer law to depict the flow pattern. They found that the porosity parameter increases the magnitude of the temperature and concentration of particles. Seddeek^[9] probed convective heat transfer in a fluid immersed in a porous medium. Analytical results for the nonlinear extended Darcy flow were described by Jha and Kaurangini^[10]. Bakar et al.^[11] studied the convective heat mode in the stagnant flow with an improved Forchheimer relation. Aziz et al.^[12] computed the traveling wave solution for time-dependent viscoelastic fluid by way of a flat absorbent plate.

Non-Newtonian fluids are highly advantageous fluids because of their intrinsic flow features and ability to overcome the difficulties encountered by various field experts at different stages. Non-Newtonian fluids are relied upon extensively in various industrial processes because of their adjustable factors, and include biological fluids (such as synovial, blood, and salvia), clinical formulations, cosmetics and toiletries, paints, food items (such as jellies, jams, and soups), and multi-phase combinations (such as gas-liquid mixtures, oil-water emulsions, butter, and froths and foams). It is notable that the extensive utility and applicability of non-Newtonian fluids are due to the nonlinear expression between the applied shear stress and the resulting strain. This flexible and most complex behavior of non-Newtonian fluids allows it to work in and benefit numerous processes. In order to understand these complexly structured liquids, they are categorized into three main classes: shear thinning (pseudoplastic), shear thickening (dilatant), and thixotropic fluids. It is clear from the literature that various fluid models have been proposed to explicate the transport features of these three classes of fluids. However, in the current work, the Sutterby fluid model is proposed, which is dual in nature and exposes the characteristics of both pseudoplastic and dilatant materials. Guha and Pradhan^[13] explored the comportment of convectional power-law fluid flow by planar configuration. Bijjan et al.^[14] demonstrated the outcome of thermal transportation on a streamed flow of dilatant liquid in a cylinder. Azhar et al.^[15] surveyed the stagnant Sutterby nanofluid flow with entropy scrutiny. Hayat et al.^[16] investigated the aspects of peristaltic Sutterby fluid flow with the Roseland heat flux in a vertical confinement. In another study, Hayat et al.^[17] revealed the peristaltic flow in a curved channel with Joule heating and radiative heat flux. Xie and Jian^[18] investigated the power-law MHD fluid flow through a micro-parallel spinning channel. Thermophysical scrutiny of three-dimensional MHD nanoparticles over an exponential stretched surface was analytically studied by Rehman et al.^[19]. Dual solutions for the stagnation point flow of Carreau fluid were numerically determined by Hashim et al.^[20]. The flow assisting and opposing impacts were scrutinized by using stability analysis. Hayat et al.^[21] disclosed the peristaltic study of Sutterby liquid flow through porous media. In another study, Hayat et al.^[22] demonstrated the impression of wall properties on the peristaltic flow of the radiative Sutterby fluid in curved surfaces. The stretched flow of Burger's fluid was analyzed by Soomro et al.^[23]. They computed the skin friction in the limiting case for comparative study and found excellent agreement with the published work.

The boundary-layer region in continuum mechanics is the space contained by the fluid molecules between the no-slip zone (the interaction area between the solid boundary and fluid molecules) and the free stream (the upper surface of fluid particles). The behavioral analysis of fluids within this region is of extreme essence for all researchers and fluid experts. However, in recent years, several analyses have revisited and modified the continuum fluid theory, and proposed the idea that it is possible to configure the flow phenomenon by generating slip between the boundary and fluid interaction points. These thought-provoking dynamists have also investigated whether this new prediction can help in reducing the stress drag factor on fluid molecules and improve the efficiency in motion of fluid structure within the boundary layered region. The flow problem under the slip effect in various configurations and physical circumstances is important for pervasive technological and industrial applications, especially in micro-electromechanical systems, such as the functioning of micro pumps, micro mixing of physiological samples, biological transportation, and drug delivery. Among the previously referenced studies, many researchers^[24-33] have disclosed the features of fluid phenomena by implementing the no-slip condition. However, as mentioned in Ref. [34], when the fluid is particulate, such as emulsions, suspensions, and polymeric solutions, the no-slip condition is inappropriate. In this circumstance, the sufficient condition is partial slip. Wang^[35] analyzed the flow properties over a planar stretched sheet under the application of partial slip. Anderson^[36] estimated the flow behavior of viscous fluid induced by a stretchable plate with partial slip. Matthews and Hill^[37] analyzed the impact of partial slip on the boundary layer region by constructing a comparison with the Navier velocity and linear slip conditions. Sajid et al.^[38] computed the shear stress rate affected by slip conditions for an incompressible viscous flow induced by a linearly stretched surface. A study on the partial slip flow with porosity effects over a linearly stretched sheet was conducted by Wang^[39].

MHD property is the combination of unknown magnetic properties of electrically conducting fluids. Alfeven^[40] pioneered the MHD field. The effect of magnetic fields on fluids has totally diverted the direction and attention of researchers. Thus, the analysis of fluid flow transmission with magnetic aspects has extensive utility in various engineering and industrial disciplines. These include glass blowing, metal extrusion, spin melting, heat exchanger design, fiber production, fabric sheet manufacturing, crystal growing, nuclear reactor cooling, sensors, electromagnetic casting, and magnetic drug targeting. Because of these industrial and engineering applications of magnetic fields, extensive studies have been published in this regard. Liao^[41] computed the skin friction affected by magnetic fields on the fluid flow over a linearly stretched surface. The power-law fluid in the presence of magnetic field was explored by Cortell^[42]. Recently, Usman et al.^[43] utilized a least-square procedure to analyze the transmission of heat in a channel. A decreasing attribute in the velocity profile of Newtonian fluid due to the implication of a magnetic field was shown by Ishak et al.^[44]. Haq et al.^[45] scrutinized the heat transfer mechanism of CuO-water filled in a partially heated rhombus with a heated hurdle by incorporating the FEM. In their exploration, they observed the significant influence of the inner heated square cylinder. Ellahi et al.^[46] examined the generalized Couette flow of Powell-Eyring fluid with MHD and slip effects, and found a numerical solution with a pseudo-spectral collocation scheme. Usman et al.^[47] utilized the collocation method to study the effect of thermal and velocity slip on the Casson nanofluid flow over a stretching cylinder. They considered the transverse magnetic field and considered Buongiorno's model for nanofluid analysis. In another study, Usman et al.^[48] computed the wavelet solutions to the MHD viscous fluid flow under slips and radiation effects in a channel. Their outcomes consisted of a comparative study with error and convergence analysis, and demonstrated the efficiency of the proposed algorithm. Phenomena of natural convection flow in a partially heated trapezoidal cavity loaded with a nanofluid in the presence of single-wall carbon nanotubes (SWCNTs) were analyzed using the FEM by Haq et al.^[49].

The heat transfer mechanism is one of the most important phenomena of fluid flow. Heat transfer in the system occurs because of temperature differences. This thermal change causes the heat to flow from one place to another or from the system to its surroundings. Initially, the Fourier law of heat conduction was leveraged to measure the thermal attributes of a system. Later, researchers criticized the restricted limit of Fourier's law due to its parabolic nature. They showed that the parabolic nature of Fourier's law restricts it to satisfy the causality principle. After this contradiction, researchers proposed many recommendations to overcome the flaw in Fourier's law. At last, Cattaneo^[50] was successful in his endeavor and modified Fourier's law by incorporating an additive thermal relaxation term in Fourier's equation. This amendment causes the energy equation to be hyperbolic, which allows the propagation of thermal waves within the system with a finite speed. This kind of energy transport has staggering practical applications that span from nanofluid flows to the modeling of skin burn injury by Tibullo and Zampoli^[51]. Christov^[52] transformed Cattaneo's work by replacing the time derivative with the upper convected derivative. This amendment made this law more realistic and appropriate. From then, this model has been recognized as the Cattaneo-Christov law. The structural stability and uniqueness of the Cattaneo-Christov equations were anticipated by Straughan^[53]. Han et al.^[54] analyzed heat transfer with slip conditions for the Cattaneo-Christov model. They validated their findings by constructing a comparison with a numerical scheme. Mustafa^[55] analyzed heat transfer disclosed by the Cattaneo-Christov model in the rotating flow of Maxwell fluid. Haq et al.^[56] analyzed the mixed convective boundary layer flow of Newtonian fluid induced by a stretching sheet. Hayat et al.^[57] computed the optimal solution of Darcy-Forchheimer flow of mixed-convection Maxwell nanomaterial over a heated surface. The numerical solution of radiative swirling flow over a heated rotating stretched disc by engaging homogeneous-heterogeneous reactions was explored by Khan et al.^[58]. Mahdy^[59] analyzed the bearing of variable temperature on the unsteady Casson liquid over a rotating sphere. Javaherdeh and Najjarnezami^[60] scrutinized the flow in a porous cavity by the lattice Boltzman procedure. Guo et al.^[61] explored the asymptotic solution for the asymmetric flow in a porous channel under the application of magnetic effects.

The present study is conducted to contemplate the Darcy-Forchheimer slip flow of a Sutterby fluid beside a linearly stretching sheet under magnetic and radiative heat fluxes. In light of the abovementioned literature, no one so far has attempted to elucidate this nature of the physical phenomena for the Sutterby fluid. First, the equations are modeled in the form of partial expressions by employing a boundary layer approach. Then, these constructed partial differential equations (PDEs) are converted into ordinary differential equations (ODEs) by exploiting the similarity approach and then solved computationally by the shooting method. The effects of the introduced parameters on velocity, temperature, shear stress factor, and convection rate profiles are shown realistically. Finally, comparison of the present work with previous findings is also addressed.

2 Problem structure 2.1 Model description

Consider the steady, two-dimensional electro-magnetohydrodynamic (E-MHD) boundary layer Darcy-Forchheimer slip flow of incompressible Sutterby fluid in the presence of Lorentz force^[31], which is described as

(1)

The flow is considered in a sheet that is stretched with the velocity U_w(x)=a₁^* x. A constant temperature T_w^* is anticipated at the surface of the sheet. A constant magnetic field of strength B_o is imposed in the vertical direction to the flow (see Fig. 1). We assume that the Reynolds number of the fluid is insignificant so that the magnetic field and Hall effects are negligible. We analyze the magnetic field and electric field effects on the momentum and thermal boundary layer equations. Mechanisms of heat transmission are considered by using the generalized Fourier law, nonlinear radiative heat flux, and viscous dissipation.

With these assumptions, the predominant boundary layer equations that describe the physical phenomena are defined as

(2)

(3)

(4)

(5)

Here, u₁ and v₁ are the velocity components in the x- and y-directions, respectively, M is the flow deportment index, b² is the consistency index, T^* is the temperature, ρ^* is the density, σ₁ is the magnetic permeability, E_o is the electric field strength, C_F is the Forchheimer quantity, k^* is the permeability of the absorbent medium, c_p is the particular heat, K^* is the thermal conductivity, λ₁^* is the thermal relaxation time, and is the radiative heat flux.

The significant boundary conditions that designate the physical phenomena are

(6)

Here, U_w=a₁^*x is the stretching velocity with a₁^* >0, T_w^* is the wall temperature, T_∞^* is the ambient temperature, N₀ is the slip length, and K_o represents the thermal jump length. Further, it is notable that for N₀=K_o=0, the present situation reduces to the no-slip flow.

Assume the following similarity transmutations:

(7)

Substituting Eq. (7) in Eqs. (3) and (4), we obtain

(8)

(9)

with the transformed boundary conditions of Eq. (6),

(10)

In overhead expressions, γ =λ₁^*a^* is the thermal relaxation parameter, is the magnetic parameter, is the porous medium parameter, is the electric parameter, is the inertia parameter, is the Prandtl number, is the Reynolds number, is the Deborah number, is the Eckert number, and is the radiation parameter. Mathematically, the relations for the skin friction coefficient and the local Nusselt number in the dimensionless form are expressed as

(11)

3 Numerical simulation

The resulting PDEs (3) and (4) of the present work are transformed into ODEs after applying a suitable transformation. However, these governing ODEs (8) and (9) are highly nonlinear, and it is difficult to control their analytical solutions. Keeping in mind the end goal to analyze this system with associated boundary conditions (10), we must solve it numerically. There are many numerical methods that can solve this problem, but we solve this problem by the shooting technique with the Runge-Kutta-Fehlberg (RKF) method, specifically, RKF45, because the RKF method is a technique that solves only the initial value problem. Hence, first, the governing boundary value problem is transmuted into an initial value problem. Thus, the governing equations are rewritten as

(12)

(13)

Because it is clear that the nonlinear momentum of Eq. (12) is of third order and the nonlinear energy Eq. (13) is of second order, the total order of all equations is seven. To find the computational solution to these equations by the Runge-Kutta (RK) scheme^[56], these equations should be transmuted into a system of five equations with five missing variables, by letting

(14)

under the new variables defined in Eq. (12), given by

(15)

(16)

(17)

(18)

(19)

The conforming boundary restrictions from Eq. (6) in the new constraints are defined in Eq. (10), and take the following form:

(20)

To elucidate the system of first-order differential equations presented in Eqs. (15)-(19), five initial conditions must be known, but the initial conditions at infinity are not prescribed. However, the boundary conditions are prescribed at infinity. Thus, these boundary conditions are used to generate two unknown initial conditions. Now, by letting the missing initial conditions convert the given boundary into initial constraints, new conditions are defined as

(21)

Initially, the values are chosen as -1 and 1. Now, to resolve this system of seven first-order ODEs with initial conditions, the RKF method is used. The RK fourth- and fifth-order formulae derived by Fehlberg are

(22)

where subscripts 5 and 4 denote the fifth- and fourth-order formulae, respectively. In addition, K₁, K_i, F, and Z are defined as

(23)

where i =2, 3, ..., 6.

The determination of the step size h is highly important. If we take h too large, the truncation error may be unacceptable. If h is too small, the iterative process is long enough. Thus, initially, we take the value of h as 0.1, and this is modified at each step. The coefficients in the RKF formula are defined in Table 1.

Equations (15)-(19) are solved with the fifth-order formula. The fourth-order formula is used only to estimate the truncation error, defined as

(24)

The computed values are a function of

(25)

The correct values of e yield the boundary conditions as η→∞ that satisfy the relation

(26)

where the differences between the computed and given boundary values are called residuals. If the boundary residuals are less than the error tolerance, i.e., 10^-6, it is the final solution. Moreover, Eq. (26) can be solved by using the Newton method to refine the value. This procedure is continued until it satisfies the convergence criteria (see Fig. 2).

4 Physical description

The purpose of this section is to review the properties of numerous significant parameters, such as the magnetic parameter Ha, the thermal relaxation parameter γ, the Prandtl number Pr, the thermal contamination parameter Nr, the sponginess parameter λ, the Deborah number De, the porous media parameter k₁, the Reynolds number Re, the thermal omission parameter β, the Eckert number Ec, and the momentum omission parameter δ, on the velocity f^'(η) and temperature θ (η) profiles. For this determination, Figs. 3-17 are shown.

The effect of the Deborah number De on f'(η) by fixing the other parameters is exposed through Fig. 3, because De is used in rheology to characterize the fluidity of materials under particular flow circumstances. Here, M=-0.5 is the shear thinning case, and M=0.5 corresponds to the shear thickening case. Figure 3 shows that increasing the Deborah number De accelerates the fluid velocity f'(η) in the case of shear thickening circumstances, and an opposite behavior is observed for the shear diminishing case. The behavior of the sponginess parameter k₁ on f^'(η) is shown in Fig. 4 in the presence (Ha≠ 0) and absence of the magnetic parameter (Ha=0). Increasing the values of the sponginess parameter k₁ decelerates the fluid velocity f^'(η). Physically, the resistance of the porous medium increases through permeability, which causes a reduction in the fluid velocity. It also shows that in the presence of Ha, there is a noticeable reduction in the fluid velocity as compared with the absence of Ha. Figure 5 shows the influence of the velocity slip δ on the fluid velocity f^'(η). It is determined that the fluid velocity decreases with the increasing slip factor. From Fig. 5, it is likely that the dimensionless velocity f^'(η) is restrained under no-slip condition (δ =0). Figure 6 shows the behavior of the magnetic parameter Ha on f^'(η) in the presence of velocity slip (δ ≠ 0) and in the absence of velocity slip (δ=0). As the magnetic parameter Ha increases, the magnetic force normal to the flow direction dominates the viscous influence. Consequently, the velocity in the axial direction decreases with the increasing magnetic parameter Ha. It also shows that a reduction in the velocity occurs more rapidly with slip. Figure 7 shows the effects of the Reynolds number Re on the velocity profile for M=0.5 and M=-0.5 because the Reynolds number is an important dimensionless quantity in fluid mechanics used to describe flow configurations in diverse fluid flow circumstances. Greater values of the Reynolds number relate to an increase in inertial forces as compared with viscous forces, which increase the fluid velocity f^'(η) for the shear thickening case M=0.5 as compared with the shear thinning case M=-0.5. The impacts of the Eckert number Ec with and without the velocity slip δ on f^'(η) are shown in Fig. 8.

The Eckert number is a dimensionless value used in flow calculations. It represents the relation between the kinetic energy and enthalpy of the flow, and is used to illustrate degeneracy. The comportment of the skin friction coefficient versus the velocity slip parameter δ by varying the Hartmann number Ha for classical and Sutterby fluids is shown in Fig. 9. It is revealed that the skin friction is higher for the Sutterby fluid case, i.e., for M=De=Re=0.5, as compared with the classical case, i.e., for M=De=Re=0. It is notable that the present flow for the Sutterby fluid model reduces to the Newtonian fluid flow by considering M=De=Re=0. In addition, it shows that because of the increase in Ha, the skin friction decreases. Figure 10 shows the behavior of the skin friction versus the Hartmann number Ha by considering the shear thickening (M>0) and shear thinning (M < 0) cases for Re=2 and Re=3. It is clear from this figure that for increasing M and Re, the skin friction is reduced.

The effects of the thermal radiation parameter Nr in the absence (β =0) and in the presence of the thermal slip (β ≠0) on the temperature profile θ (η) are shown in Fig. 11. It is revealed that the radiation parameter causes an increase in the fluid temperature, and the associated layer thickness increases as well. From the figure, it is clear that for β =0, the temperature solution satisfies the boundary condition, as well as in the presence of thermal slip, by keeping in mind the concept of heat transfer coefficient. The behaviors of thermal slip (β ≠ 0) for Pr =0.5 and Pr =1.0 are shown in Fig. 12 by fixing the other parameters. It is shown that the thermal slip friction in fluid particles decreases, which causes a reduction in the fluid temperature. The influence of the Eckert number Ec on the dimensionless temperature θ (η) in the absence (β =0) and presence of thermal slip (β ≠ 0) is shown in Fig. 13. Increasing the value of the Eckert number Ec increases the temperature θ (η). The absence of thermal slip (β =0) shows that the temperature θ (η) is higher as compared with the presence of thermal slip (β ≠ 0). Physically, this means that an increase in the external magnetic field causes an increase in the kinetic energy of the fluid particles. This increase in the kinetic energy entails an increase in Joule heating (ohmic dissipation) and viscous degeneracy. Subsequently, the fluid temperature increases.

Figure 14 shows the influence of thermal slip, velocity slip, and Prandtl number on the heat transfer coefficient. It shows that the heat transfer coefficient increases for higher values of the Prandtl number, and it decreases with the velocity and thermal slip. The impact of the thermal relaxation parameter γ on the fluid temperature is shown in Fig. 15. Higher values of γ correspond to a decrease in heat due to the reduction in the fluid temperature.

The impact of the thermal relaxation parameter γ in the absence of thermal slip (β =0) on the heat transfer coefficient is shown in Fig. 16 by a bar diagram.

It is observed that heat transfer increases by increasing γ. Figure 17 shows the three-dimensional portrayal of velocity.

Some physical quantities of interest are presented by varying k₁, M, Re, De, Pr, and β with the help of Tables 2 and 3. The impact of the skin friction is addressed in Table 2. It is shown that increases for k₁=0.1, 0.2, 0.3, Re=2, 3, 4, and De=0.2, 0.4, 0.6, and it decreases for M=0.1, 0.2, 0.3. The effects of different influential parameters, such as M, Pr, and β, on the heat transfer coefficient are listed in Table 3. It is observed that -θ'(0) increases by increasing Pr and M, and the opposite relation between -θ'(0) and β is shown.

5 Closing remarks

The two-dimensional E-MHD boundary layer flow of steady, constant-density Sutterby fluid by a stretching sheet with slip conditions is investigated in this study. The properties of numerous flow adjusting parameters on the dimensionless velocity and temperature are scrutinized. From the present inquiries, the conclusions can be summarized as follows.

(ⅰ) The Hartmann number condenses the non-dimensional velocity profile f^'(η) and increases the fluid temperature.

(ⅱ) Increasing the value of the velocity slip parameter δ corresponds to a decrease in the non-dimensional velocity profile f^'(η).

(ⅲ) Increasing the value of the magnetic parameter Ha and the velocity slip parameter δ causes a similar comportment to the non-dimensional velocity profile f^'(η).

(ⅳ) Assisting and opposing impacts of De and k₁ are observed for increasing the values of these parameters on f^'(η).

(ⅴ) A similar impact of Re, Ha, δ, k₁, and Ec is observed for the fluid velocity.

(ⅵ) Increasing the values of the thermal relaxation parameter γ and thermal slip β reduces the fluid temperature, whereas an opposite behavior of the radiation parameter Nr is observed.

Acknowledgements  The authors would like to express their gratitude to King Khalid University, Abha 61413, Saudi Arabia for providing administrative and technical support.