EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5463 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Simulation for Obtaining Advanced Engineering Applications Through an Important Fluid Problem Ali H. Tedjani1,∗, Mohamed M. Khader1 1 Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia Abstract. The current work examines the transfer of heat/mass in a thin liquid film (TLF) with high viscosity on a surface that is continuously stretching. Through thorough mathematical mod- eling, a system of PDEs is established. The application of appropriate similarity transformations results in the formulation of a set of nonlinear ODEs from the original system of PDEs. Here, we introduce a reliable numerical technique to analyze the solution characteristics of the proposed problem. This technique depends on applying shifted Chebyshev polynomials of the sixth-order (SCP6s). Initially, we obtain an approximate expression for the derivatives by using SCP6s. In this study, we employ a spectral collocation procedure, leveraging the unique characteristics of SCP6s. This approach enables the transformation of the presented problem into a set of nonlinear algebraic equations. Notably, we emphasize the presentation of specific theorems that pertain to examining the convergence of the proposed approximation technique. Our present findings exhibit a remarkable concurrence with the investigations that have been previously reported. 2020 Mathematics Subject Classifications: 41A30, 65N12, 76F12 Key Words and Phrases: Heat and mass fluxes, Casson fluid, Liquid thin film, Porous medium, Heat generation, Chebyshev-spectral collocation technique ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5463 Email addresses: Ahtedjani@imamu.edu.sa (A. H. Tedjani), mmkhader@imamu.edu.sa (M. M. Khader) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 2 of 18 Nomenclature Symbol Description Symbol Description a, b Positive constants B0 Magnetic field strength C Fluid concentration Cr Reference concentration Cfx Drag force coefficient d Positive constant D Brownian diffusion coefficient E Activation energy parameter Ea Coefficient of activation energy Ec Eckert number h(t) Film thickness k Porous medium permeability K1 Chemical reaction coefficient K2 Chemical reaction parameter M Magnetic field parameter Nux Nusselt number Pr Prandtl number qc Mass flux coefficient qt Heat flux coefficient Q Heat generation coefficient S Unsteadiness parameter Sc Schmidt number T Temperature Tr Reference temperature us Sheet velocity x, y Cartesian coordinates ρ Fluid density β Casson parameter µ Fluid viscosity σ Electrical conductivity coefficient κ Thermal conductivity coefficient η Similarity variable ν0 Kinematic viscosity at the slit δ Heat generation parameter θ Non-dimensional temperature λ Porous parameter ψ The dimensional stream function α Viscosity parameter ϵ Slip velocity coefficient ϕ Non-dimensional concentration ε1 Thermal conductivity parameter ε2 Diffusivity parameter ε3 Temperature difference parameter ε4 Slip velocity parameter γ Dimensionless film thickness 1. Introduction The phenomenon of flow obtained through a stretching surface has been a subject of great interest and complexity, captivating the attention of engineers and researchers for many years. This phenomenon finds applications in various industrial sectors, en- compassing processes such as coating, painting, melt-spinning, and metal production. In these operations, a continuous stretching of a thin sheet or film of material occurs, in- ducing flow in the initially stationary fluid. This flow behavior is intricate and relies on the characteristics of both the material being stretched and the fluid itself. Gaining a thorough understanding of this phenomenon is crucial for optimizing and enhancing the efficiency of these processes [10]. Further, enhancing efficiency and performance further requires a grasp of heat transfer mechanisms and their applications in various engineering systems. As such, prior research [22] provides a full explanation of this phenomenon and its important applications. Film flow refers to the movement of a slender liquid layer over a surface, driven by external forces such as gravity or shear stress. It is characterized by the presence of at A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 3 of 18 least one free boundary. Thin films are employed in the field of coating. Coating plays a crucial role in numerous manufacturing processes for microelectronic devices, including hard discs, and microelectronic chips. The model took into account various factors, in- cluding inertial forces, disjoining friction, surface tension, and curvature, that influence the system dynamics. The significance of studying thin film flow (TFF) over an unsta- ble stretching sheet (SS) has motivated numerous researchers ([5], [8]) to conduct studies addressing various practical problems within this field. The analysis of varying heat flux and its effects on Casson thin film flow and HT generated by an unsteady SS is highly significant. This inquiry gives a pivotal role in overseeing the complexities of heat mass transfer, particularly in contexts associated with cooling processes. Megahed’s research [18], which underscores the detailed examination of the influence of variable heat flux, offers valuable insights into the intricacies of Casson TFF and the dynamics of HT associated with it. Moreover, the significance of fluid properties exhibiting variability and their influence on Casson thin film behavior can be elucidated through the research conducted by Mahmoud and Megahed [16]. Their study delves into the intricate aspects of how fluid properties, subject to variations, play a crucial role in shaping the characteristics of Casson thin films. This research contributes to the broader comprehension of fluid dynamics, shedding light on the interplay of various factors influencing the motion and thermal properties of Casson liquid films. It is important to highlight that analytical solutions for most ODEs can not be ob- tained. As a result, many researchers strive to obtain approximate solutions for ODEs. Here, we specifically concentrate on numerically solving the proposed model. With this goal in mind, a hybrid spectral collocation method (SCM) is applied, employing the SCP6s. It is important to emphasize that SCM has been effectively utilized to computationally handle various mathematical models. In our proposed technique, the utilization of these basis functions results in a system of nonlinear algebraic equations. The solution of these equations is used to expand the functions in a polynomial form. The primary merit of this methods, as opposed to previously developed SCMs in references ([7], [11], [12], [14]), is its distinctive feature. Another benefit lies in the convergence order, which indicates the method’s higher accuracy compared to others. Additionally, we will give an error analysis of the Chebyshev approximations used in this technique. From the literature mentioned earlier, it has been observed that the combined influence of heat generation, viscous dissipation, magnetohydrodynamic mass HT chemical reaction, and variable fluid properties in a porous medium has not been explored in the context of a Casson TLF flow across an unsteady stretching surface. Hence, the motivation and novel aspect of this study lies in utilizing the SCM, specifically depending on SCP6s to investigate the problem under study. 2. Flow analysis This work inspects the behavior of a liquid film flowing over a stretched sheet immersed in a porous medium, taking into consideration the endurance of a magnetic field. The flow in the thin film, with a thickness denoted h(t), is induced by the stretching of the sheet A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 4 of 18 along the x-axis at a velocity of us(x, t). The precise definition of this phenomenon can be found in the reference cited as [15]. us(x, t) = b x(1− at)−1, a, b > 0. It must seen that for the current analysis, only quantities of t that satisfy t < 1 a are considered acceptable. In addition, we consider the incompressible flow within a TLF that is situated within a porous medium with permeability k (see Figure 1 for more detail). Figure 1. Geometry of flow problem Moreover, in our investigation, we assume that the sheet is related to variable heat flux qt(x, t) and variable mass flux qc(x, t). Both of these fluxes depend on x & t and are related according to the equation provided in reference [15]: qt(x, t) = −κ ( ∂T ∂y ) = x2 d Tr (1− at) −5 2 , (1) qc(x, t) = −D ( ∂C ∂y ) = x2 dCr (1− at) −5 2 . (2) In the previous equations, the construction of equations (1) and (2) recognizes that both the heat flux and the mass flux between the SS and the LF increase proportionally with x2. Furthermore, these fluxes exhibit a time-dependent growth pattern, indicating that both the heat flux and mass flux intensify over time. In our analysis, we deem the re- spectable interrelationship between both µ(T ) and κ(T ) of the Casson fluid and the fluid temperature. This connection is captured by the formula presented in reference [19]: µ(T ) = µ0 Exp [ −α ( T − T0 Tr d x2 (1−at)2 1 κ0 √ ν0 b )] , κ(T ) = κ0 [ 1 + ε1 ( T − T0 Tr d x2 (1−at)2 1 κ0 √ ν0 b )] , D(C) = D0 [ 1 + ε2 ( C − C0 Cr d x2 (1−at)2 1 D0 √ ν0 b )] . A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 5 of 18 The parameters in these equations play a decisive role in characterizing the attitude of the model. In this work, our focus is to maintain a similarity solution. To achieve this objective, the transverse applied magnetic fieldB(t) is expressed by the approach presented by Abel et al. [3]: B(t) = B0 [1− at] −1 2 , where B0 is a constant. Considering the previously stated assumptions, the temperature, concentration, and velocity, fields in the thin liquid layer are believed by the momentum, concentration, and thermal energy equations in 2Dim., which will be formulated as follows ([1], [19]): ∂u ∂x + ∂v ∂y = 0, (3) ∂u ∂t + u ∂u ∂x + v ∂u ∂y = 1 ρ ( 1 + 1 β ) ∂ ∂y ( µ(T ) ∂u ∂y ) − µ(T ) ρk ( 1 + 1 β ) u− σB2 0 ρ u, (4) ∂T ∂t + u ∂T ∂x + v ∂T ∂y = 1 ρcp ∂ ∂y ( κ(T ) ∂T ∂y ) + Q ρcp + µ(T ) ρcp ( 1 + 1 β )( ∂u ∂y )2 , (5) ∂C ∂t + u ∂C ∂x + v ∂C ∂y = ∂ ∂y ( D(C) ∂C ∂y ) −K1(C − C0). (6) Here, we must note that in Thumma and Mishra’s previously published study [24], the last term of the energy equation indicates viscous dissipation, which explains the conversion of kinetic energy into thermal energy due to viscous forces inside the fluid. Also, the expression for the term Q is derived from the work of Gomathy and Kumar [9]: Q = κus xν0 δ(T − T0). In the fluid model, the term λ is > 0 when heat is generated and < 0 when heat is absorbed. The boundary conditions (B.Cs) suitable for the current system are defined as follows [19]: v = 0, u = us(x, t), −κ∂T ∂y = qt(x, t), −D∂C ∂y = qc(x, t), at y = 0, (7) v = dh dt , ∂C ∂y = 0, ∂u ∂y = 0, ∂T ∂y = 0, at y = h(t). (8) We will define the newly introduced dimensionless variables f , ϕ, and θ, along with the similarity parameter η, in the following manner [19]: ψ = √ ν0b (1− at) −1 2 xf(η), η = √ b ν0 (1− at) −1 2 y, (9) ϕ(η) = C − C0 Cr d x2 (1−at)2 1 D0 √ ν0 b , θ(η) = T − T0 Tr d x2 (1−at)2 1 κ0 √ ν0 b . A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 6 of 18 In this context, the variable ψ represents the physical stream function, which inherently satisfies the mass conservation equation (3). Moreover, u & v are dependent on ψ, as illustrated by the subsequent equations: u = ∂ψ ∂y , v = −∂ψ ∂x . The mathematical problem presented in the system (4)-(6) along with B.Cs (7)-(8) can now be simplified to a set of ODEs and their corresponding B.Cs, which are given below:( 1 + 1 β )( f ′′′ − α θ′f ′′ ) e−α θ + ff ′′− f ′2−Mf ′−S (η 2 f ′′ + f ′ ) −λ ( 1 + 1 β ) f ′ e−α θ = 0, (10) 1 Pr ( (1 + ε1 θ)θ ′′ + ε1 θ ′2)+ fθ ′ − 2f ′ θ − S ( 2θ + η 2 θ ′ ) + δ Pr (1 + ε1θ) θ + Ec ( 1 + 1 β ) f ′′ 2e−α θ = 0, (11) (1 + ε2 ϕ)ϕ ′′ + ε2 ϕ ′2 + Sc [ fϕ ′ − 2f ′ϕ− S ( 2ϕ+ η 2 ϕ′ ) −K2ϕ ] = 0, (12) f(0) = 0, f ′(0) = 1, θ′(0) = −1 1 + ε1θ(0) , ϕ′(0) = −1 1 + ε2ϕ(0) , (13) f(γ) = γ 2 S, f ′′(γ) = 0, θ′(γ) = 0, ϕ′(γ) = 0. (14) The explanations for each parameter that impacts the momentum, concentration, and energy fields are outlined in the following form: S = a b , M = σB2 0 ρ b2 , λ = ν0(1− at) kb , Ec = κ0 b 5 2 cp d Tr √ ν0 , Pr = µ0 cp κ0 , Sc = ν0 D0 , K2 = K1(1− at) b . Furthermore, equation (9) provides a means to compute the parameter γ, given by the following expression: γ = ( bρ µ0(1− at) ) 1 2 h(t). In the complete system of equations (10)-(12), it is evident that γ (unknown parameter) needs to be determined. 3. Industrial quantities There are three considerable physical measures that play a crucial role in practical applications. These quantities must be acknowledged and discussed individually. They A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 7 of 18 include Nux, Cfx, and Shx. The definitions and expressions for these quantities are provided below: CfxRe 1 2 = − ( 1 + 1 β ) e−α θ(0)f ′′(0), NuxRe −1 2 = 1 θ(0) , ShxRe −1 2 = −ϕ′(0), where Re = us(x,t)x ν0 is the local Reynolds number. 4. Approximation functions and their implementation 4.1. Some properties of the SCP6s and its approximations We introduce some of the basic definitions and properties of the shifted Chebyshev polynomials of the sixth kind (SCP6s) ([6], [13], [23]) to suit their use in solving the proposed model in the interval [0, ℏ]. The shifted CP6s on [0, ℏ], may be constructed by using the linear transformation z = (2/ℏ)η − 1 as Tk(η) = Tk((2/ℏ)η − 1) [17]. We can get the polynomials {Tk(η)}∞k=0 with the help of the following recurrence relation [2]: Tk(η) = ((2/ℏ)η − 1)Tk−1(η)− ( (−1)k(2k + 1) + k(k + 1) + 1 4k(k + 1) ) Tk−2(η), k = 2, 3, ... , where T0(η) = 1, T1(η) = (2/ℏ)η − 1. The analytic form of Tk(η) is defined by: Tk(η) = k∑ i=0 ci,k η i, (15) where ci,k = 22i−k ℏi(2i+ 1)!  ∑⌊ k 2 ⌋ r=⌊ i+1 2 ⌋ (−1)0.5k+r+i(2r+i+1)! (2r−i)! , if k even,∑⌊ k−1 2 ⌋ r=⌊ i 2 ⌋ 2(−1)0.5(k+1)+r+i(r+1)(2r+i+2)! (k+1)(2r−i+1)! , if k odd. The function ψ(η) ∈ L2[0, ℏ] will be approximated as a finite series sum with the first (m+ 1)-terms in the following form: ψm(η) = m∑ ℓ=0 ζℓ Tℓ (η). (16) Theorem 1. [4] The achieved error εm = |ψ(η)−ψm(η)| in approximating the function ψ(η) by ψm(η) which is defined in (16) will be estimated with the criterion: εm ≤ 2−m. A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 8 of 18 Theorem 2. Let ψ(η) be estimated by CP6s as (16), then: ψ(n) m (η) ≃ m∑ i=n i∑ k=n ζi χi, k,n η k−n, n > 0, (17) where χi, k,n = (k!/(k − n)!) ci,k, but ci,k is defined in (15). 4.2. Implementation the SCM-SCP6s We use the SCP6s to estimate the solution of the given problem (10)-(12), a set of basis functions with the SCM. We implement SCM to solve numerically the model (10)-(12). We estimate f(η), θ(η), and ϕ(η) by fm(η), θm(η), and ϕm(η), respectively as follows: fm(η) = m∑ i=0 ai Ti(η), θm(η) = m∑ i=0 bi Ti(η), ϕm(η) = m∑ i=0 ci Ti(η). (18) By collecting the Eqs.(10)-(12), with the formulae (17)-(18), we find:( 1 + β−1 ) ( f (3)m (η)− α θ(1)m (η) f (2)m (η) ) e−α θm(η) + fm(η) f (2)m (η)− ( f (1)m (η) )2 −M ( f (1)m (η) ) − S ( 0.5 η f (2)m (η) + f (1)m (η) ) − λ ( 1 + β−1 ) ( f (1)m (η) ) e−α θm(η) = 0, (19) 1 Pr ( (1 + ε1 θm(η)) θ(2)m (η) + ε1 ( θ(1)m (η) )2 ) + fm(η) θ(1)m (η)− 2f (1)m (η) θm(η) − S ( 2θm(η) + 0.5η θ(1)m (η) ) + δ Pr (1 + ε1θm(η)) θm(η) + Ec ( 1 + β−1 ) ( f (2)m (η) ) e−α θm(η) = 0, (20) (1 + ε2 ϕm(η))ϕ(2)m (η) + ε2 ( ϕ(1)m (η) )2 + Sc [ fm(η)ϕ(1)m (η)− 2 f (1)m (η)ϕm(η) − S ( 2ϕm(η) + 0.5 η ϕ(1)m (η) ) −K2 ϕm(η) ] = 0. (21) These equations (19)-(21) can be collocated with m−2 of points ηp (the roots of Tm−2(η)) as follows:( 1 + β−1 ) ( f (3)m (ηp)− α θ(1)m (ηp) f (2) m (ηp) ) e−α θm(ηp) + fm(ηp) f (2) m (ηp)− ( f (1)m (ηp) )2 −M ( f (1)m (ηp) ) − S ( 0.5 η f (2)m (ηp) + f (1)m (ηp) ) − λ ( 1 + β−1 ) ( f (1)m (ηp) ) e−α θm(ηp) = 0, (22) A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 9 of 18 1 Pr ( (1 + ε1 θm(ηp)) θ (2) m (ηp) + ε1 ( θ(1)m (ηp) )2 ) + fm(ηp) θ (1) m (ηp)− 2f (1)m (ηp) θm(ηp) − S ( 2θm(ηp) + 0.5η θ(1)m (ηp) ) + δ Pr (1 + ε1θm(ηp)) θm(ηp) + Ec ( 1 + β−1 ) ( f (2)m (ηp) ) e−α θm(ηp) = 0, (23) (1 + ε2 ϕm(ηp))ϕ (2) m (ηp) + ε2 ( ϕ(1)m (ηp) )2 + Sc [ fm(ηp)ϕ (1) m (ηp)− 2 f (1)m (ηp)ϕm(ηp) − S ( 2ϕm(ηp) + 0.5 η ϕ(1)m (ηp) ) −K2 ϕm(ηp) ] = 0. (24) Also, the B.Cs (13)-(14) may be expressed by substituting from Eq.(18) in (13)-(14) to get the following equations: m∑ i=0 ai = 0, ( m∑ i=0 bi θ (1) m (0) )( 1 + ε1 m∑ i=0 (−1)i bi ) = −1, m∑ i=0 ai T (1) i (0) = 1, ( m∑ i=0 ci T (1) i (0) )( 1 + ε2 m∑ i=0 (−1)i ci ) = −1, (25) m∑ i=0 ai Ti(γ) = 0.5γ S, m∑ i=0 ai T (2) i (γ) = 0, m∑ i=0 bi T (1) i (γ) = 0, m∑ i=0 ci T (1) i (γ) = 0. (26) Eqs.(22)-(26), assign a system of 3(m+1) algebraic equations. We use the Newton iteration technique to solve these equations for the unknowns ai, bi, ci, i = 0, 1, ...,m. Hence we compute the residual error function (REF) for the approximation formulas (17) and (18) in the following forms [21]: REF[f(η)] = ( 1 + β−1 ) ( f (3)m (η)− α θ(1)m (η) f (2)m (η) ) e−α θm(η) + fm(η) f (2)m (η)− ( f (1)m (η) )2 −M ( f (1)m (η) ) − S ( 0.5 η f (2)m (η) + f (1)m (η) ) − λ ( 1 + β−1 ) ( f (1)m (η) ) e−α θm(η), REF[θ(η)] = 1 Pr ( (1 + ε1 θm(η)) θ(2)m (η) + ε1 ( θ(1)m (η) )2 ) + fm(η) θ(1)m (η)− 2f (1)m (η) θm(η) − S ( 2θm(η) + 0.5η θ(1)m (η) ) + δ Pr (1 + ε1θm(η)) θm(η) + Ec ( 1 + β−1 ) ( f (2)m (η) ) e−α θm(η), REF[ϕ(η)] = (1 + ε2 ϕm(η))ϕ(2)m (η) + ε2 ( ϕ(1)m (η) )2 + Sc [ fm(η)ϕ(1)m (η)− 2 f (1)m (η)ϕm(η) − S ( 2ϕm(η) + 0.5 η ϕ(1)m (η) ) −K2 ϕm(η) ] . A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 10 of 18 5. Validation of the numerical solution To test the precision of the numerical results obtained through the SCM-SCP6s method, we present a table of data that offers a comparative evaluation between our current find- ings and the previously published results by Pal and Saha [20]. This tabular data (Table 1) allows for a detailed analysis and comparison of the accuracy between the two sets of results solutions. The acquired solutions show a strong degree of consistency with the formerly published data. It is clear that the latest results closely match the former data published, demonstrating a strong correlation and validating our methodology’s accuracy. We used the Mathematica 11 to write and remedy all the codes on Dell Inspiron 15 (3593) Workstation, Processor: 11th Gen Intel(R) Core(TM) i7-1165G7 and 2.80GHz 1.69 GHz, 1 TB storage, and 32 GB Ram DDR3. Table 1. Comparison of −f ′′ (0) and γ for various values of S using SCM-SCP6s method when M = λ = α = 0 and β → ∞ SSSpt¡ - Spt¿ Pal and Saha [20] Current work γ −f ′′ (0) γ −f ′′ (0) 1.2 1.4 1.6 1.8 1.1277820 0.8210349 0.5761762 0.3563941 1.442627 1.012785 0.642402 0.309146 1.12778195 0.82103481 0.57617616 0.35639402 1.44262680 1.01278475 0.64240092 0.30914589 6. Interpretation of numerical results In this work, we investigated the dynamic behavior of a TLF flowing over a solid sur- face that could be stretched in two dimensions, using the Casson model. A porous medium was used to contain this surface, and as a result, the fluid’s characteristics were variable. Furthermore, we studied the effect of the heat generation, magnetic field, chemical reac- tion, and viscous dissipation on the system. We analyzed how these factors influenced the velocity, heat distribution, and mass transfer rate. To visualize the variations in the calcu- lated results, we employed the SCM-SCP6s method and utilized Mathematica software to generate graphical representations. Figure 2 presents f ′(η), θ(η), and ϕ(η) corresponding to different unsteadiness parameters S. The figure demonstrates a consistent pattern: as S increase, there is a consistent upward trend observed in the free stream concentration ϕ(γ), free stream velocity f ′(γ), and fluid velocity f ′(η). However, there is a contrasting trend observed in θ(η), sheet temperature θ(0), and film thickness γ, which decrease as the unsteadiness parameters increase. Physically, a reduction in the unsteadiness parameter suggests that the flow becomes more stable and less erratic over time. This diminished unsteadiness can contribute to a slower and more prolonged spreading of the nanofluid on the surface, resulting in an augmentation of the thin film thickness. A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 11 of 18 0.5 1.0 1.5 2.0 Η 0.2 0.4 0.6 0.8 1.0 f 'HΗL S=1.5, Γ=0.8643 S=1.2, Γ=1.3322 S=0.8, Γ=2.3187 Β=0.5, Λ=0.5 M=0.5, Α=0.1 0.5 1.0 1.5 2.0 Η 0.2 0.4 0.6 S=1.5 S=1.2 S=0.8 ∆=0.2, Pr=3.0 Ec=0.2, K2=0.1 ¶1=0.2, ¶2=0.2 ΘHΗL ΦHΗL S=0.8, 1.2, 1.5 Figure 2. (a) f ′(η) for various S (b) ϕ(η) and θ(η) for various S The significance of λ on θ(η), ϕ(η), and f ′(η) is clearly depicted in Figure 3. Through- out the entire range of η values, the presence of λ leads a noticeable reduction in both the film thickness γ, and f ′(η). This reduction occurs because the porous medium hinders the flow. The ϕ(η) and θ(η) fields exhibit an opposite trend compared to the film thickness and fluid velocity when considering the same porous parameter λ. This behavior leads to an enhancement in both the free surface concentration ϕ(γ) and the free surface tempera- ture θ(γ). Physically, a decrease in λ signifies a diminished porous quality of the medium, indicating a more permissive environment for fluid flow. Consequently, the nanofluid en- counters reduced resistance while moving through the porous medium. This decreased hindrance promotes a more effective dispersion of the nanofluid on the surface, resulting in an augmentation of the thin film thickness. 0.5 1.0 1.5 Η 0.2 0.4 0.6 0.8 1.0 f 'HΗL Λ=1.0, Γ=1.1276 Λ=0.5, Γ=1.3322 Λ=0.0, Γ=1.7140 Β=0.5, S=1.2 M=0.5, Α=0.1 0.5 1.0 1.5 2.0 Η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ∆=0.2, Pr=3.0 Ec=0.2, K2=0.1 ¶1=0.2, ¶2=0.2 ΘHΗL ΦHΗL Λ=1.0, 0.5, 0.0 Figure 3. (a) f ′(η) for various λ (b) ϕ(η) and θ(η) for various λ Figure 4 presents the profiles of f ′(η), ϕ(η), and θ(η) for various quantities of M . It is important to note that the involvement of M gives rise to the Lorentz force. This force acts as a resistance to the flow, leading to a decrease in both f ′(η), and γ. Physically, a A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 12 of 18 reduction in the magnetic parameter points a weaker magnetic field affect on the nanofluid flow. Consequently, the magnetic forces restraining the flow become less influential. This diminished magnetic effect results in decreased resistance, leading to a thicker thin film. Simultaneously, the diminishing magnetic parameter relaxes the magnetic forces, causing an increase in the nanofluid velocity. The weakened magnetic field allows for freer fluid flow, enhancing the overall velocity of the thin film. 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Η 0.2 0.4 0.6 0.8 1.0 f 'HΗL M=1.0, Γ=1.2515 M=0.5, Γ=1.3322 M=0.0, Γ=1.43065 Β=0.5, S=1.2 Λ=0.5, Α=0.1 0.5 1.0 1.5 2.0 Η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ∆=0.2, Pr=3.0 Ec=0.2, K2=0.1 ¶1=0.2, ¶2=0.2 ΘHΗL ΦHΗL M=0.0, 0.5, 1.0 Figure 4. (a) f ′(η) for various M (b) ϕ(η) and θ(η) for various M Figure 5 gives the influence of the Casson parameter β on f ′(η), ϕ(η), and θ(η) dis- tributions. The graph illustrates that as β increases, both ϕ(γ) and ϕ(0) show a rise, while a reverse trend occurs for both θ(0) and θ(η). This indicates an enhancement in the distribution of ϕ(η) throughout the film layer. Additionally, it is seen that both f ′(η) and γ grow as β gets closer to zero. Physically, β represents the yield stress in the Casson fluid model. A decrease in the Casson parameter implies a lower yield stress, indicating reduced resistance to flow initiation. This allows the nanofluid to spread more easily over the surface, leading to a thicker thin film. Simultaneously, the decrease in the Casson pa- rameter signifies a shift to a more fluid-like behavior, reducing resistance within the thin film and resulting in an increased velocity of the nanofluid as it encounters less opposition to its motion. 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Η 0.2 0.4 0.6 0.8 1.0 f 'HΗL Β=2.0, Γ=1.0509 Β=1.0, Γ=1.1670 Β=0.3, Γ=1.47685 M=0.5, S=1.2 Λ=0.5, Α=0.1 0.5 1.0 1.5 2.0 Η 0.2 0.4 0.6 ∆=0.2, Pr=3.0 Ec=0.2, K2=0.1 ¶1=0.2, ¶2=0.2 ΘHΗL ΦHΗL Β=0.3, 1.0, 2.0 Figure 5. (a) f ′(η) for various β (b) ϕ(η) and θ(η) for various β A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 13 of 18 Figure 6 shows how ε1 and δ parameters affect the modes of HT. A material’s capacity to conduct heat is measured by its thermal conductivity. It measures how quickly heat may move through a material. Heat will be absorbed and conducted more efficiently by a fluid with a high thermal conductivity. This implies that when heat is generated inside or transferred to a fluid, the fluid will swiftly disperse the heat throughout its volume. The fluid temperature will consequently rise quickly as obvious from the figure for high thermal conductivity values. Also, the heat generation parameter δ describes how much heat is generated inside the fluid per unit volume or unit time. Increases in the heat generation parameter indicate that the fluid is producing more heat. The fluid’s temperature rises as a result of the heat produced as indicated in Figure 6(b), which increases the fluid’s overall energy content. The temperature elevation is rooted in the fact that the extra heat produced within the nanofluid enhances the thermal environment, intensifying the overall thermal conditions. This increased thermal effect causes the temperature distribution to expand throughout the thin film, resulting in elevated temperatures across the surface. 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.05 0.10 0.15 0.20 0.25 0.30 0.35 ΘHΗL ∆=0.2, Pr=3.0 Ec=0.2, K2=0.1 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, ¶2=0.2 ¶1=0.0, 0.5, 1.0 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.1 0.2 0.3 ΘHΗL ¶1=0.2, Pr=3.0 Ec=0.2, K2=0.1 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, ¶2=0.2 ∆=0.0, 0.5, 1.0 Figure 6. (a) θ(η) for various ε1 (b) θ(η) for various δ Figure 7 illustrates the variations in temperature profiles θ(η) with respect to changes in Ec and Pr. Clearly, θ(η), θ(0), and θ(γ) decrease as Pr grows because Pr represents the ratio of the fluid’s momentum diffusivity to its thermal conductivity. Additionally, the graph demonstrates that as Ec increases, θ(η), θ(0), and θ(γ) increase, while the di- mensionless thickness of the thin film remains unchanged. Physically, the rationale for the augmented temperature distribution can be traced to the heightened efficiency in convert- ing kinetic energy into heat associated with the elevated Eckert number. This increased efficiency in converting energy results in a more pronounced thermal effect within the Casson nanofluid. Consequently, this intensified thermal impact prompts the expansion of θ(η) throughout the thin film, leading to elevated temperatures across the surface. A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 14 of 18 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ΘHΗL ∆=0.2, Pr=3.0 ¶1=0.2, K2=0.1 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, ¶2=0.2 Ec=0.0, 0.5, 1.0 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.1 0.2 0.3 0.4 0.5 ΘHΗL ∆=0.2, Ec=0.2 ¶1=0.2, K2=0.1 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, ¶2=0.2 Pr=1.0, 3.0, 6.2 Figure 7. (a) θ(η) for various Ec (b) θ(η) for various Pr Figure 8 depicts the concentration profile ϕ(η) as it varies with changes in both the chemical reaction parameter K2 and the diffusion parameter ε2. According to observa- tions, the concentration field shrinks as the diffusion parameter rises close to the sheet, but the opposite tendency is seen farther away from the sheet, and the film thickness stays the same. In the presence of K2, the dominant trend in ϕ(η), ϕ(γ), and ϕ(0) is a diminish- ing behavior. Physically, the underlying physical cause for the diminishing concentration distribution can be attributed to the increased chemical reaction parameter, which encour- ages more substantial chemical reactions. This heightened reactivity leads to a depletion of the nanofluid concentration, observed uniformly across the thin film. Consequently, the concentration distribution experiences a decrease as a result of this reduction in nanofluid concentration. 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ΦHΗL ∆=0.2, Pr=3.0 ¶1=0.2, K2=0.1 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, Ec=0.2 ¶2=0.0, 0.5, 1.0 0.2 0.4 0.6 0.8 1.0 1.2 Η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ΦHΗL ∆=0.2, Ec=0.2 ¶1=0.2, Pr=3.0 M=0.5, S=1.2 Λ=0.5, Α=0.1 Β=0.5, ¶2=0.2 K2=0.0, 0.2, 0.4 Figure 8. (a) ϕ(η) for various ε2 (b) ϕ(η) for various K2 A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 15 of 18 Table 2. Values of Re 1 2 xCfx, Re −1 2 x Nux and Re −1 2 x Shx for assorted values of S, λ,M, β, ε1, ε2, δ, Ec, Pr and K2 with α = 0.1 and Sc = 0.5 S λ M β ε1 δ Ec Pr ε2 K2 Re 1 2 xCfx Re −1 2 x Nux Re −1 2 x Shx 0.8 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.18835 2.35897 1.31384 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.09541 2.75015 1.40679 1.5 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 2.76277 3.16960 1.47903 1.2 0.0 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 2.45820 2.99581 1.47490 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.09541 2.75015 1.40679 1.2 1.0 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.61594 2.56201 1.33933 1.2 0.5 0.0 0.5 0.2 0.2 0.2 3.0 0.2 0.1 2.89612 2.82487 1.42999 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.09541 2.75015 1.40679 1.2 0.5 1.0 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.28230 2.68161 1.38368 1.2 0.5 0.5 0.3 0.2 0.2 0.2 3.0 0.2 0.1 4.01456 2.61875 1.43951 1.2 0.5 0.5 1.0 0.2 0.2 0.2 3.0 0.2 0.1 2.36545 2.85957 1.35455 1.2 0.5 0.5 2.0 0.2 0.2 0.2 3.0 0.2 0.1 1.97450 2.91741 1.30505 1.2 0.5 0.5 0.5 0.0 0.2 0.2 3.0 0.2 0.1 3.09521 2.69494 1.40676 1.2 0.5 0.5 0.5 0.5 0.2 0.2 3.0 0.2 0.1 3.09573 2.82625 1.40683 1.2 0.5 0.5 0.5 1.0 0.2 0.2 3.0 0.2 0.1 3.09619 2.93935 1.40689 1.2 0.5 0.5 0.5 0.2 0.0 0.2 3.0 0.2 0.1 3.09588 2.77670 1.40680 1.2 0.5 0.5 0.5 0.2 0.5 0.2 3.0 0.2 0.1 3.09451 2.70993 1.40673 1.2 0.5 0.5 0.5 0.2 1.0 0.2 3.0 0.2 0.1 3.09279 2.64161 1.40660 1.2 0.5 0.5 0.5 0.2 0.2 0.0 3.0 0.2 0.1 3.10994 3.61736 1.40781 1.2 0.5 0.5 0.5 0.2 0.2 0.5 3.0 0.2 0.1 3.07391 2.02823 1.40527 1.2 0.5 0.5 0.5 0.2 0.2 1.0 3.0 0.2 0.1 3.03892 1.41776 1.40276 1.2 0.5 0.5 0.5 0.2 0.2 0.2 1.0 0.2 0.1 3.05675 1.75811 1.40443 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.1 3.09541 2.75015 1.40161 1.2 0.5 0.5 0.5 0.2 0.2 0.2 6.2 0.2 0.1 3.10717 3.44181 1.40061 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.0 0.1 3.09541 2.75015 1.35394 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.5 0.1 3.09541 2.75015 1.47372 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 1.0 0.1 3.09541 2.75015 1.56344 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.0 3.09541 2.75015 1.40202 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.2 3.09541 2.75015 1.42802 1.2 0.5 0.5 0.5 0.2 0.2 0.2 3.0 0.2 0.4 3.09541 2.75015 1.46936 In Table 2, the drag force, HT coefficient, and mass transfer coefficient are expressed in non-dimensional form as CfxRe 1 2 x , NuxRe −1 2 x , and ShxRe −1 2 x , respectively. These pa- rameters are calculated against different quantities of the controlling parameters. It can be easily observed that the non-uniformity of the heat flux has a notable impact on the quantities of the Nusselt number. Specifically, increasing the porous parameter, magnetic number, heat generation parameter, and Ec tends to decrease the Nusselt number. Con- versely, the reverse tendency is shown for all other parameters, with the exception of the A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 16 of 18 diffusion and chemical reaction parameters, leading to a boost in SFC. Furthermore, the observed decrease in the Sherwood number values for increasing magnetic number, Casson parameter, thermal conductivity parameter, heat generation parameter, and Eckert num- ber, is accompanied by an increasing trend for the Sherwood number with the unsteadiness parameter, diffusion parameter, and chemical reaction parameter. These trends can be attributed to the non-uniformity of the mass flux applied in this study. To further confirm the efficiency, accuracy and effectiveness of the proposed method, the REF values were calculated across different values of the order of approximation m under the same boundary conditions and the same parameter values mentioned above and are listed in Table 3. Also, we can improve the accuracy by increasing m. Table 3. The REF of the solution via various m. η m REF − f(η) REF − θ(η) REF − ϕ(η) 4 5.6512E− 05 3.0591E− 05 5.8593E− 05 0.0 7 3.0256E− 08 5.2553E− 07 3.5873E− 08 4 1.1082E− 05 9.5244E− 05 1.2564E− 05 0.2 7 0.0831E− 08 0.2953E− 08 2.9208E− 07 4 3.1197E− 05 3.1167E− 04 4.5573E− 05 0.4 7 8.3404E− 08 6.2091E− 07 6.2205E− 08 4 2.0063E− 05 3.5503E− 05 8.0764E− 05 0.6 7 1.0412E− 08 8.5497E− 07 7.2183E− 08 4 0.0907E− 05 1.5569E− 05 4.5273E− 05 0.8 7 6.1523E− 07 4.5708E− 08 1.2894E− 08 4 9.2297E− 05 8.2130E− 05 0.1046E− 05 1.0 7 1.9923E− 08 3.5708E− 08 5.0076E− 07 7. Conclusions This study is focused on the analysis of the proposed important problem. The SCM- SCP6s technique is utilized to solve a system of ODEs that govern the proposed model. Different physical parameters’ effects on the computed findings for velocity, heat, and mass transfer rate are examined and represented through tables and figures. The velocity profile explains a diminishing bearing with increasing magnetic number and Casson number. While the free surface quantities show an increase, the film thickness falls as the Casson parameter, magnetic number, and unsteadiness increase. As the Eckert parameter and heat generation parameter are increased, the sheet temperature and temperature profiles rise, which raises the free surface temperature. On the other hand, Pr exhibits the opposite pattern, where an increase causes the sheet temperature, temperature profiles, and free surface temperature to decrease. The rates of heat and mass transport are accelerated by the unsteadiness and thermal conductivity parameters. Future research in this field could explore hybrid non-Newtonian nanofluid thin film flow under varying density, mass and heat fluxes, and Ohmic heating. A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 17 of 18 References [1] Z. Abbas, T. Hayat, M. Sajid, and S. Asghar. Unsteady flow of a second-grade fluid film over an unsteady stretching sheet. Mathematical and Computer Modelling, 48:518–526, 2008. [2] W. M. Abd-Elhameed and Y. H. Youssri. Sixth-kind chebyshev spectral approach for solving fractional differential equations. Int. J. Nonlinear Sci. Numer. Simul., 20:191–203, 2019. [3] M. S. Abel, J. Tawade, and M. M. Nandeppanavar. Effect of non-uniform heat source on mhd heat transfer in a liquid film over an unsteady stretching sheet. International Journal of Non-Linear Mechanics, 44:990–998, 2009. [4] A. G. Atta, W. M. Abd-Elhameed, G. M. Moatimid, and Y. H. Youssri. Advanced shifted sixth-kind chebyshev tau approach for solving linear one-dimensional hyper- bolic telegraph-type problems. Math. Sci., 6:1–15, 2022. [5] R. C. Aziz, I. Hashim, and S. Abbasbandy. Flow and heat transfer in a nanofluid thin film over an unsteady stretching sheet. Sains Malaysiana, 47:1599–1605, 2018. [6] J. P. Boyd. Chebyshev and Fourier Spectral Methods. Dover, New York, USA, 2nd edition, 2000. [7] M. Delkhosh and H. Cheraghian. An efficient hybrid method to solve nonlinear differential equations in applied sciences. Comp. Appl. Math., 41:1–15, 2022. [8] A. Elham and A. M. Megahed. Mhd dissipative casson nanofluid liquid film flow due to an unsteady stretching sheet with radiation influence and slip velocity phenomenon. Nanotechnology Reviews, 11:463–472, 2022. [9] G. Gomathy and B. R. Kumar. Variable thermal conductivity and viscosity effects on thin film flow over an unsteady porous stretching sheet. Journal of Porous Media, 14:77–94, 2023. [10] P. Gupta and A. Gupta. Heat and mass transfer on a stretching sheet with suction or blowing. The Canadian Journal of Chemical Engineering, 55:744–746, 1977. [11] M. M. Khader. Numerical study for unsteady casson fluid flow with heat flux using a spectral collocation method. Indian Journal of Physics, 96:777–786, 2021. [12] M. M. Khader, A. Eid, and M. Adel. Implementing the vieta-lucas collocation op- timization method for mhd casson and williamson’s model under the effects of heat generation and viscous dissipation. Journal of Mathematics, 2022:1–13, 2022. [13] M. M. Khader and K. M. Saad. On the numerical evaluation for studying the frac- tional kdv, kdv-burger’s, and burger’s equations. The European Physical Journal Plus, 133:1–13, 2018. [14] M. M. Khader and R. P. Sharma. Evaluating the unsteady mhd micropolar fluid flow past stretching/shirking sheet with heat source and thermal radiation: Implementing fourth order predictor-corrector fdm. Math. Comput. Simulat., 181:333–350, 2021. [15] I. C. Liu, A. M. Megahed, and H. W. Hung. Heat transfer in a liquid film due to an unsteady stretching surface with variable heat flux. ASME Journal of Applied Mechanics, 80:12–23, 2013. [16] M. A. A. Mahmoud and A. M. Megahed. Mhd flow and heat transfer characteristics in A. H. Tedjani, M. M. Khader / Eur. J. Pure Appl. Math, 18 (1) (2025), 5463 18 of 18 a casson liquid film towards an unsteady stretching sheet with temperature-dependent thermal conductivity. Brazilian Journal of Physics, 47:512–523, 2017. [17] J. C. Mason and D. C. Handscomb. Chebyshev Polynomials. Chapman and Hall, CRC, New York, NY, Boca Raton, 2003. [18] A. M. Megahed. Effect of slip velocity on casson thin film flow and heat transfer due to an unsteady stretching sheet in the presence of variable heat flux and viscous dissipation. Applied Mathematics and Mechanics, 36:1273–1284, 2015. [19] I. G. Nourhan and A. M. Megahed. Hydromagnetic nanofluid film flow over a stretch- ing sheet with prescribed heat flux and viscous dissipation. Fluid Dynamics and Material Processing, 18:1373–1388, 2022. [20] D. Pal and P. Saha. Analysis of unsteady magnetohydrodynamic radiative thin liquid film flow, heat and mass transfer over a stretching sheet with variable viscosity and thermal conductivity. International Journal for Computational Methods in Engineer- ing Science and Mechanics, 22:400–409, 2021. [21] K. Parand and M. Delkhosh. Operational matrices to solve nonlinear volterra- fredholm integrodifferential equations of multi-arbitrary order. Gazi University Jour- nal of Science, 294:895–907, 2016. [22] M. D. Shamshuddin, S. R. Mishra, O. A. Beg, and A. Kadir. Adomian decomposition method simulation of von Kármán swirling bioconvection nanofluid flow. Journal of Central South University, 26:2797–2813, 2019. [23] M. A. Snyder. Chebyshev Methods in Numerical Approximation. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1966. [24] T. Thumma and S. R. Mishra. Effect of viscous dissipation and joule heating on magnetohydrodynamic jeffery nanofluid flow with and without multi-slip boundary conditions. Journal of Nanofluids, 73:516–526, 2018.