Fourier heat conduction model is the basic key for thermoelectric transport models [1]. However, a fundamental drawback of the widely used Fourier's law is that it describes a parabolic energy equation, which is valid for instantaneous propagation of thermal disturbances and consequently, the energy transport is anisotropic and assessment of unrealistics infinite heat propagation speed. To address this deficiency, a general structure known as Fourier's law with thermal relaxation was developed by Cattaneo [2] where a thermal relaxation time, representing the finite thermal conduction speed, was assigned. Exploiting this idea, Christov [3] goes a step further and improves the model by framing it in a frame-indifferent way using Oldroyd's theory of upper-convected derivatives, so that compatibility with the tenets of rational thermodynamics is achieved. The earlier work of Christov [3] initiated a series of studies of heat transfer with non-Fourier effects. Straughan [4] deduced the variation in thermal results subect to natural convective flow by improsing the non-tradtional Fourier's approach.

At the same time, there is a great interest in investigating different types of non-Newtonian fluids, ranging from everyday products to industrial liquids such as blood, ketchup, cheese, drilling mud, and colloidal suspensions. Of a variety of constitutive models, the 1944 Eyring–Powell fluid, originally proposed by Powell and Eyring, has been recognized for its ability to replicate Newtonian-like behavior at both low and high shear rates. Its empirical nature provides a more realistic description of complex rheological behavior, particularly for shear thickening and thinning phenomena. Several studies, including those by [5, 6, 7], have validated the versatility of the Eyring–Powell model in diverse flow configurations. [8] further advanced this field by incorporating ferromagnetic effects into the Eyring–Powell model while utilizing the Cattaneo-Christov (CC) heat flux theory. Their work highlights the necessity of coupling advanced heat conduction models with non-Newtonian flow theory to achieve accurate predictions of coupled heat and momentum transport. Recent progress in non-Newtonian fluid mechanics has expanded these insights: [9] investigated viscoelastic nanofluid transport under temperature gradients, while [10] analyzed Casson polymeric flows with radiative and non-Fourier effects. For bioconvective applications, [11] examined Eyring–Powell nanofluids over Riga surfaces, and [12] explored hybrid nanofluids with Hall current effects. Additionally, [13] developed numerical methods for fractional advection-diffusion systems, and [14, 15] addressed chemically reactive flows and double-diffusion phenomena in micropolar nanofluids.

Recent developments in fluid dynamics have emphasized the importance of varying sheet thickness on the behaviour of non-Newtonian fluids especially in industrial processes involving calendering and coating. Zahid et al. [16] studied the calendering of non-isothermal visco-elastic materials and reported that temperature dependent viscosity has a major influence on the variable sheet thickness and final properties of the material in the process. Similarly, Anwar et al. [17] studied non-Newtonian flow of fluids over stretching sheets in porous space and stressed on the effect of thermal conductivity and magnetic fields in changing the pattern of flow and heat transfer rates. Abbas et al. [18] realized the curved surface flow with chemical specieis for non-Newtonian fluids. Nonlinear stretching surfaces play an important impact in transport phenomenon, especially for complex materials. It is reported through some investigative studies that magnetic field, buoyancy and convective boundary conditions play a vital role in flow features and heat transfer characteristics in such systems [19, 20, 21, 22]. These observations are also useful for the industrial processes regarding the nonlinear stretching surface. The Modified Eyring–Powell model is particularly useful in replicating shear-thinning behaviour without developing singularities at zero shear rate as seen in models such as power-law or Bingham. Modified Eyring–Powell model has recently been utilized to consider many engineering problems, such as temperature dependent viscosity in web coating process, non-similar boundary layer flow and heat transfer over a wedge, and chemically reacting nanofluid flow over wedges. These applications demonstrate the model's practical applicability and accuracy in describing the rheological response of the non-Newtonian fluid in various engineering and biomedical practices (see [23, 24, 25]).

Motivated by the above scenario, the present research is pursued to address the impacts of non-Fourier heat flux in Darcy–Forchheimer flow of a Modified Eyring–Powell fluid towards a nonlinear stretching surface with variable thickness. Thermal results are claimed under the assumptions of temperature dependent thermal conductivity. The problem is simplified into dimensionless forms. The homotopy analysis scheme, with excellent accuracy is followed to entertain the simulations.

We investigate a two-dimensional Darcy–Forchheimer flow of a modified Eyring–Powell fluid on a non-linear stretching sheet with variable thickness. The stretching surface starts with a slit placed at the origin with x-axis placed along the sheet and y-axis perpendicular to it. The thermal conductivity is considered as temperature-dependent, and heat diffusion is described using a non-Fourier (Cattaneo–Christov) heat flux model. The surface profile is y=A⁢(x+b)1−n2. The boundary layer assumptions are introduced and the equations of continuity, momentum, and energy written out in the boundary layer coordinate system.

∂ x u + ∂ y v = 0

u ⁢ ∂ x u + v ⁢ ∂ y u = μ ρ ⁢ ∂ y ⁢ y u + q ρ ⁢ ( 1 ρ ⁢ B ⁢ C ) q ⁢ ( 1 + q + 2 6 ⁢ c 2 ⁢ ( ∂ y u ) 2 ) × ( ∂ y u ) q − 1 ⁢ ∂ y ⁢ y u − μ ρ ⁢ k ∗ ⁢ u − C b k ∗ ⁢ u 2

u ∂ x T + v ∂ y T + λ { ( u ( ∂ x u ) ( ∂ x T ) + v ( ∂ y u ) ( ∂ y T ) ) + v ∂ y u ( ∂ x T ) + 2 u v ∂ x ⁢ y T + u 2 ∂ x ⁢ x T + v 2 ∂ y ⁢ y T } = − 1 ρ ⁢ c p ⁢ k ⁢ ( T ) ⁢ ∂ y ⁢ y T

with conditions

u = U w ⁢ ( x ) = U 0 ⁢ ( x + b ) n , v = 0 , T = T w ⁢ at ⁢ y = A ⁢ ( x + b ) 1 − n 2 , u → 0 , T → T ∞ ⁢ as ⁢ y → ∞

The thermal conductivity is assumed to be k⁢(T)=k∞⁢(1+ϵ⁢Θ), where k∞ is the thermal conductivity of the ambient fluid, ϵ symbolizes a small coefficient for thermal conductivity, and Θ is dimensionless temperature. The constants (b, c, d) have no dimension, and μ is the dynamic viscosity and ρ the fluid density. The factors B, C, and q are related to the Modified Eyring–Powell (MEP) fluid model. Where cp is the specific heat, T is the fluid temperature, T∞ the ambient temperature, and u,v are the velocity components in the x, y direction, respectively. Following quantities are adopted for simplification:

χ = 2 n + 1 ⁢ U 0 ⁢ ν ⁢ ( x + b ) n + 1 ⁢ F ⁢ ( ξ ) , ξ = n + 1 2 ⁢ ( U 0 ν ) ⁢ ( x + b ) n − 1 , u = U 0 ⁢ ( x + b ) n ⁢ F ′ ⁢ ( ξ ) , v = − n + 1 2 ⁢ U 0 ⁢ ν ⁢ ( x + b ) n − 1 ⁢ ( F ⁢ ( ξ ) + ξ ⁢ n − 1 n − 1 ⁢ F ′ ⁢ ( ξ ) ) , Θ ⁢ ( ξ ) = T − T ∞ T w − T ∞

By the use of Eq. (5) into the Eqs.(2-4), we get

( 1 + K ⁢ q ⁢ n + 1 2 ⁢ δ ⁢ F ′′ ) q − 1 ⁢ ( 1 − q + 2 3 ⁢ δ ⁢ ( F ′′ ) 2 ) ⁢ F ′′′ + F ⁢ F ′′ − 2 ⁢ n n + 1 ⁢ ( F ′ ) 2 − 2 n + 1 ⁢ β ⁢ F ′ − 2 n + 1 ⁢ f ⁢ r ⁢ ( F ′ ) 2 = 0

( 1 + ϵ Θ ) Θ ′′ + ϵ Θ ′ ⁣ 2 + P r ( F Θ ′ + γ ( n − 3 2 F Θ ′ F ′ − n + 1 2 F 2 Θ ′′ ) ) = 0

with

F ⁢ ( 0 ) = 1 , F ⁢ ( 0 ) = 1 − n 1 + n , Θ ⁢ ( 0 ) = 1 , F ′ ⁢ ( ∞ ) = 0 , Θ ⁢ ( ∞ ) = 0

Now we have F⁢(ξ)=f⁢(ξ−α)=f⁢(η) and Θ⁢(ξ)=θ⁢(ξ−α)=θ⁢(η), then Eq. (7-9) becomes

( 1 + K ⁢ q ⁢ n + 1 2 ⁢ δ ⁢ f ′′ ) q − 1 ⁢ ( 1 − q + 2 3 ⁢ δ ⁢ ( f ′′ ) 2 ) ⁢ f ′′′ + f ⁢ f ′′ − 2 ⁢ n n + 1 ⁢ ( f ′′ ) 2 − 2 n + 1 ⁢ β ⁢ f ′ − 2 n + 1 ⁢ f ⁢ r ⁢ ( f ′ ) 2 = 0

( 1 + ϵ Θ ) Θ ′′ + ϵ Θ ′ ⁣ 2 + P r ( f Θ ′ + γ ( n − 3 2 f Θ ′ f ′ − n + 1 2 f 2 Θ ′′ ) ) = 0

with

f ⁢ ( 0 ) = 1 , f ⁢ ( 0 ) = 1 − n 1 + n , Θ ⁢ ( 0 ) = 1 , f ′ ⁢ ( ∞ ) = 0 , Θ ⁢ ( ∞ ) = 0

Here K(=1μ⁢B⁢C) and δ(=U03⁢(x+b)3⁢n−12⁢ν⁢c2) illustrating the fluid parameters, f⁢r(=2⁢cb⁢(x+b)k∗⁢(n+1)),β(=2⁢νk∗⁢(n+1)⁢U0⁢(x+b)n−1),γ the thermal relaxation parameter and P⁢r(=να). The wall drag force is:

C f = 2 ⁢ τ w ρ ⁢ U w 2

In which

τ w = μ ⁢ ∂ y u + 1 B ⁢ C ⁢ ∂ y u − 1 6 ⁢ B ⁢ C 3 ⁢ ( ∂ y u ) 3 at y = A ⁢ ( x + b ) 1 − n 2

The dimensionless form is

R ⁢ e x ⁢ C f 2 = n + 1 2 ⁢ ( 1 + K − n + 1 2 ⁢ K ⁢ δ 3 ⁢ f ′′ ⁢ ( 0 ) ) ⁢ f ′′ ⁢ ( 0 )

where R⁢ex=Uw⁢(x+b)ν the local Reynolds number.

The convergent solutions of Eqs. (10) and (11) are solved using the HAM, in the presence of the boundary conditions (5)-(7), (12). The solutions are given in series, consistent with the nature of the problem. The initial guesses and associated linear operators are given by:

f 0 ⁢ ( η ) = 1 − e − η , θ 0 ⁢ ( η ) = e − η , L f = f ′′′ − f ′ , L θ = θ ′′ − θ ′

Such that

L f ⁢ [ C 1 + C 2 ⁢ e η + C 3 ⁢ e − η ] = 0 , L θ ⁢ [ C 4 ⁢ e η + C 5 ⁢ e − η ] = 0 .

in which Ci⁢(i=1−5) elucidate the arbitrary constants. An important part of the HAM method is the H-curve, which presents the convergence of the Eqs. (10) and (11) subject to Eqs. (12). The optimal ranges for the convergence-control parameters are −1.85≤hf≤−0.15 and −1.6≤hθ≤−0.4. The H-curve derived is shown in Figure 1. Table 1 presents accuracy of HAM simulations. The illustration of residual error for HAM has been plotted in Figure 2.

Figure 1.

H-Curve for hf and hθ.

Figure 3.

K via f′.

This work aims to provide a comprehensive analysis of the Darcy–Forchheimer flow of a Modified Eyring–Powell material. Updated heat flux approach is followed to model the problem. The noteworthy observations are:

The velocity decreases as the porosity parameter and Forchheimer number are increased.