Finite element analysis of thin film flows in non-Newtonian third grade Sisko fluid for lifting and drainage applications

Abstract

The study of thin film-like flow has emerged as a critical area in fluid dynamics, driven by its relevance to numerous applications in industries such as semiconductor manufacturing, food processing, irrigation systems, surface coating, and oil refining. This study presents a finite element method to investigate the behavior of third-grade non-Newtonian Sisko fluids in lifting and drainage scenarios. By extending the formulations of third-grade fluids, Galerkin’s finite element method is particularly developed to derive solutions and conduct a comprehensive analysis. The analysis focuses on the impact of various parameters, including film thickness, liquid velocity profile, volumetric flux, viscosity, shear stress, gravitational effects, density, and boundary conditions. The reliability, consistency, and accuracy of the results are assessed by evaluating residuals for different scenarios. The accuracy of the results obtained by the proposed framework is higher than the previous studies as the residuals of proposed Galerkin’s FEM is found as low as ${10}^{-17}$ , whereas traditional semi-analytical techniques stay around ${10}^{-6}$ to ${10}^{-2}$ . Results also reveal that the gravitational parameter inversely affects the velocity profile during lifting, while the material constant exhibits a direct relationship. Conversely, in the drainage scenario, the material constant has an inverse effect on the velocity profile. To the best of our knowledge, the thin film flow of third-grade fluid in both lifting and drainage scenarios has not been attempted using Galerkin’s FEM so far.

Introduction

Numerous natural phenomena, including lava, water-filled eyes, and raindrops on windows, can be observed as thin film flows. When a fluid follows a vertical item in a way that maintains the shape of the object and viscous forces, this phenomenon is known as free drainage. ¹ Industrial uses for these flows include oil refining procedures, paint finishing, building and public works, chip manufacturing, and laser cutting.^2–4 Based on Newtonian fluids, the first thin-film research was done in Brien and Schwartz. ³ While this process is long-lasting, it is insufficient for nonlinear analysis of non-Newton liquids, including blood, lubricants with polymeric additives, melted plastic gels, and meals like honey and ketchup.^5,6 Siddiqui et al.^7,8 used third-grade and Phan-Thein-Tanner fluids that flow along an inclined surface to alleviate the drainage issues in elevation. Siddiqui et al. ⁹ analyzed a thin-film flow involving vertical cylinders and fourth-grade fluids. Alam et al. ¹⁰ studied pseudoplastic fluid thin films. Non-Newtonian fluid flow via a circular tube was examined by Deiber and Santa Cruz. ¹¹ Yih ¹² conducted the initial investigations on laminar flows in free surfaces to flow types. Analyses of turbulent flows have been added by Landau ¹³ and Stuart. ¹⁴ Surface tension was taken into consideration when performing stability analyses by Nakaya ¹⁵ and Lin. ¹⁶ Hydrothermal investigation of a hybrid nanofluid on a vertical plate with slip effects was carried out by Zangooee et al. ¹⁷ Several characteristics of magneto-hyperbolic-tangent type liquids were examined by Gulzar et al. ¹⁸ Nanofluid in a vertical tube with a polynomial border was examined by Najafabadi et al. ¹⁹ Nayak et al. ²⁰ examined numerically mixed convection nanofluid past an isothermal thin-needle metallic nanomaterial. The significance of Lorentz forces on radiative nanofluid concerning various limitations was studied by Algehyne et al. ²¹ Boundary layer flow of non-Newtonian fluids containing planktonic microorganisms was solved by Zaher et al. ²² In a hybrid nanofluid, Sara et al. ²³ examined the thin bloodstream using electroosmotic forces. Yadav et al. ²⁴ and his co-researchers analyzed the pulsatile flow of Jeffrey fluid through an inclined overlapped stenosed artery. Ashraf et al. ²⁵ used Adomian Decomposition Method (ADM) to investigate the lifting and drainage of a Sisko fluid film over a vertically moving cylinder, considering surface tension effects.

In general, nonlinear systems governing non-Newtonian Sisko fluids lack an exact solution or are difficult to solve analytically. To investigate the solutions to such nonlinear problems, researchers have proposed approximate numerical approaches. Numerical approaches have also been proposed to solve the thin film flow of third-grade fluids with different kinds of movements, such as vertical, horizontal, and angular. These include Finite Difference Methods (FDMs),^26–30 Homotopy Perturbation Method (HPM),^31,32 Homotopy Asymptotic Method (HAM), ³³ Adomian Decomposition Method (ADM), ³⁴ Optimal Homotopy Asymptotic Method (OHAM), ³⁵ variational iteration method ³⁶ and Differential Transformation Method (DTM). ³⁷ Due to the emergence of modern sophisticated computational techniques like neural networks^38,39 and evolutionary methods,^40–43 numerical methods are now considered a basic tool.

Thin film flows of non-Newtonian Sisko fluids are observed in several natural phenomena like tear films, blood flows, lava flows and crude oil flows. They are effectively used in modern industrial scenarios, including coating processes and oil refining. Primary studies mostly addressed Newtonian fluids, but numerous real-world liquids such as polymeric solutions, blood, and food materials show non-Newtonian behavior. Therefore, more sophisticated modeling techniques are needed. Main focus of previous studies has been on investigating thin film flows of non-Newtonian Sisko fluids using approximate FDMs and perturbation based semi-analytical analytical techniques. However, these techniques are often limited by norms of simple geometries, weak nonlinearity, or convergence issues. A particular computing tool that can be used for obtaining approximate solutions to many complex problems that other techniques would find intractable is the Finite Element Method (FEM).^44,45 The comprehensive framework of FEM entails splitting the solution domain into a limited number of finite elements or simple subdomains. Despite the proven benefits of the FEM in treating vigorous nonlinearity and complex boundaries, its application to lifting and drainage scenarios of thin film flows of third-grade Sisko fluids, remains unexplored. Applications of the Galerkin’s FEM^44–46 based on weighted-residual formulation can be found in a wide range of fields, including electromagnetics, heat transport, fluid dynamics, and solid mechanics. Some particular instances of its applications include finding solutions of second order obstacle problems,^47–49 a system of third-order ⁵⁰ and fourth-order ⁵¹ obstacle problems, singular two-point boundary value problems, ⁵² and third-grade fluid flow down an inclined plane. ⁵³ There are applications of Galerkin’s FEM that arise in areas as diverse as solid mechanics, heat transfer and electromagnetism, ⁵⁴ and fluid dynamics. ⁵⁵ Finite elements have been utilized in various ways to solve boundary value problems.

In the realm of solid mechanics, it is a reputable numerical method. There are several approaches to solving boundary-value problems that make use of finite elements. In certain formulations, variational approaches are taken into consideration, while in others a weighted residual approach is used. Due to its versatility, adaptability to the geometry of the solution space, and accuracy—especially in complex engineering and scientific problems—FEM is more accurate as compared to other numerical methods. Since FEM is so good at managing complicated geometries, boundary conditions, and material heterogeneity, it is an effective tool for solving issues in the real world where irregular domains are frequent. Whereas FDMs are dependent on rigid grids and may not be able to handle a variety of forms or material qualities. FEM employs discretization into smaller, more manageable sub-domains, or elements that may be tailored to match any geometry with more accuracy and flexibility. The analytical approaches known as perturbation methods depend on sensitive parameters, which might make them less applicable to situations in which these parameters are ill-defined or nonexistent. Due to their approximation nature, these approaches may not yield precise outcomes for highly nonlinear or complicated systems unless the problem meets the presumptions inherent in the methodology. Furthermore, there is no guarantee that these approaches will always converge. Conversely, FEM is independent of minute parameters or particular presumptions on the nonlinearity of the system. It can efficiently handle very nonlinear systems by breaking the issue down into more manageable components when local approximations are performed. As a result, FEM is more reliable for a larger variety of problems. Furthermore, compared to lower-order approximations employed in other approaches, FEM’s ability to use higher-order shape functions provides greater accuracy. It is an excellent option for resolving challenging boundary value issues in science and engineering because of its solid mathematical foundations, which guarantee accurate answers.

This study attempts to address the above-mentioned research gap by developing a robust Galerkin’s FEM framework to precisely model as well as analyze thin film flows of third-grade non-Newtonian Sisko fluids with lifting and drainage scenarios. The present work seeks to answer the following particular research questions: How can Galerkin’s FEM be efficiently utilized to model complex thin film flows of third-grade and Sisko fluids? What roles do material nonlinearity, gravitational forces, and surface effects play in determining flow behavior? How do the Galerkin’s FEM solutions compare in accuracy and stability to previous approximate analytical results? What physical interpretations can be drawn regarding flow stability, and film thickness under various operating conditions?

The manuscript is organized as follows: Section “Description of thin film flow of non-Newtonian third-grade Sisko fluid” presents the mathematical formulation of non-Newtonian third-grade Sisko fluid. Section “General framework of Galerkin’s finite element” presents a general framework of Galerkin’s finite element method for solving boundary value problems. Section “Galerkin’s finite element application to non-Newtonian third-grade Sisko fluid flow” extends the solution procedure of Galerkin’s finite element method for thin film flow of a third-grade Sisko fluid by using linear Lagrange polynomials as the elements and weight functions. Section “Results and discussions” provides the solution found by the proposed FEM framework. A detail discussion and comparative analysis of the obtained results for the lifting and drainage cases are also presented in the same section. Section “Significance and novelty of the current study” describes the contributions of this study. In the end, Section “Conclusions” provides the conclusions and some future directions.

Description of thin film flow of non-Newtonian third-grade Sisko fluid

A few presumptions are made to simplify the Navier–Stokes equations because of the thin film. Among these presumptions are:

• The film thickness is significantly less than the flow’s typical length scale.

• Differences in the flow along the film’s thickness are far more important than those in the lateral direction.

• External pressure gradients and surface forces, such as surface tension, are what mainly propel the flow.

Often, these approximations lead to the lubrication approximation, which reduces the complexity of the momentum and continuity equations by assuming that pressure is almost constant throughout the film thickness and that vertical (normal to the surface) velocities are significantly smaller than horizontal velocities. Thin film flow refers to the behavior of fluid layers with small thicknesses compared to their lateral dimensions, occurring in various natural and industrial processes like liquid spreading, coating, lubrication, and biological systems like tear films in the eye. Studying thin film flow dynamics involves understanding fluid behavior under forces like gravity, surface tension, viscous stresses, and external pressure gradients. Due to the small thickness of the fluid layer, simplifying assumptions and mathematical approximations can be made, leading to more tractable models. Thin film flow dynamics are primarily governed by Navier–Stokes equations, which describe viscous fluid motion. The foundational equations governing the dynamics of thin film flow are articulated as follows^7,8:

$$\mathrm{div}\mathbf{V}=0.$$

(1)

Equation (1) expresses the divergence of the velocity field ( $\mathbf{V}$ ), ensuring the conservation of mass. For incompressible fluids, the divergence of the velocity field $\mathbf{V}$ is zero, according to this continuity equation. It guarantees that the flow of fluid maintains mass, which means that the rate at which fluid enters and exits a region is equal.

$$\rho [\frac{\partial \mathbf{V}}{\partial t}+(\mathbf{V}.\nabla )\mathbf{V}]=\nabla .\mathbf{T}+\rho \mathbf{b}.$$

(2)

Equation (2) examines the temporal evolution of the velocity field V, focusing on the interaction between inertial terms, the gradients of the Cauchy stress tensor T, and an external force density $\rho$ . It describes fluid substance motion, and balancing forces on a fluid element, including the rate of change of velocity, convective acceleration, stress tensor divergence, and body force $\mathbf{b}$ , accounting for both inertial and external effects.

$$\mathbf{T}=-p\mathbf{I}+\mathbf{S}.$$

(3)

Here, the stress tensor is $\mathbf{S}$ , the pressure is $p$ , and the unit tensor is $\mathbf{I}$ . Equation (3) describes the Cauchy stress tensor $\mathbf{T}$ as the product of pressure $p$ and the extra stress tensor $\mathbf{S}$ , where $\mathbf{S}$ represents additional stress contributions.

$$\mathbf{S}=[\mathrm{a}+\mathrm{b}{\left|\sqrt{\frac{1}{2}\mathrm{tr}({{\mathit{A}}_1}^2)}\right|}^{n-1}{\mathit{A}}_1].$$

(4)

Equation (4) expresses the extra stress tensor $\mathbf{S}$ using constants $a$ , $b$ , and $n$ , which depends on the square root of the trace of the square of the tensor ${\mathit{A}}_1$ . The material constants $a$ , $b$ , and $n$ represent the fluid’s deformation reaction parameters, while the fluid’s strain rate is represented by the square root of the trace of ${\mathit{A}}_1^2$ .

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{S}+{\mathrm{\lambda }}_1\frac{\mathrm{D}\mathbf{S}}{\mathrm{Dt}}+\frac{{\mathrm{\lambda }}_3}{2}(\mathbf{S}{\mathbf{A}}_1+{\mathbf{A}}_1\mathbf{S})+\frac{{\mathrm{\lambda }}_5}{2}\left(\mathrm{tr}\mathbf{S}\right)}\\ \multicolumn{1}{c}{{\mathbf{A}}_1=\mathrm{\mu }({\mathbf{A}}_1+{\mathrm{\lambda }}_2\frac{\mathrm{D}{\mathbf{A}}_1}{\mathrm{Dt}}+{\mathrm{\lambda }}_4{\mathbf{A}}_1^4),}\\ \multicolumn{1}{c}{{\mathbf{A}}_1=\mathbf{L}+{\mathbf{L}}^{\mathbf{T}},\mathbf{L}=\mathrm{grad}\mathbf{V}.}\end{array}$$

(5)

Equation (5) provides a constitutive connection for the extra stress tensor in complex non-Newtonian Sisko fluids, involving material constants $\mu ,{\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3\mathrm{and}{\mathrm{\lambda }}_5$ , Rivlin–Ericksen tensor, and intricate interactions. It encapsulates the nuanced mechanics of thin film flow, revealing the impact of various factors on the system’s behavior. The components (a) and (b) of Figure 1 show the geometry of lifting and drainage cases of thin film flow respectively.

Figure 1.

Geometry of (a) lifting case (b) drainage case.

Model for the lifting case

Under the assumptions of thin film lubrication, the velocity vector is initially along the $x$ -axis, therefore $\mathbf{V}=(v(x,y),0,0)$ . The pressure varies mainly with $y$ and negligibly with $x$ -axis. The dominant stress is $S_{\mathrm{xy}}$ . Due to slow flow, the convective terms are neglected. The only significant term is $\frac{d}{\mathrm{dy}}{T}_{\mathrm{xy}}$ .

The momentum equations reduce to: $x$ -direction: $\partial p/\partial x=0$ ; $y$ -direction: $\partial p/\partial y=\rho g+\partial /\partial x{T}_{\mathrm{xy}}$ . From the Sisko model for 1D simple shear flow is: $T_{\mathrm{xy}}=a\frac{\mathrm{dv}}{\mathrm{dx}}+b{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n$ . $S_{\mathrm{xy}}$ is the extra stress component, which equals $T_{\mathrm{xy}}$ in shear because pressure does not affect shear components. So, in shear components: $S_{\mathrm{xy}}=T_{\mathrm{xy}}$ . Substituting $T_{\mathrm{xy}}$ into the momentum equation and expanding the derivative we get: $-\frac{\mathrm{dp}}{\mathrm{dy}}+\rho g+a\frac{d^{2}v}{d{x}^2}+b\frac{d}{\mathrm{dx}}{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n=0$ .

By inserting the above formulations from equations (3) and (4) into equation (2), we derive:

$$\begin{array}{c}-\frac{\mathrm{dp}}{\mathrm{dx}}=0,\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{dp}}{\mathrm{dy}}+\rho g+a\frac{d^{2}v}{d{x}^2}+b\frac{d}{\mathrm{dx}}{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n=0.}\end{array}\}.$$

(6)

The momentum balancing equation (2) is used to obtain equation (6), replacing stress tensor and velocity formulations, implying a constant pressure along the $x$ -axis. This is indicated by, $-\frac{d{p}_1}{\mathrm{dx}}=0$ . Terms about gravity, such as $\rho g$ , and the second derivative of velocity, $\frac{d^{2}v}{d{x}^2}$ , are included, and the nonlinear shear stress term, which is dependent on ${\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n$ raised to the power $n$ . Further, we deduce that $p_1=p_1\left(y\right)$ ,

$$a\frac{d^{2}v}{d{x}^2}+\mathrm{nb}{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^{n-1}\frac{d^{2}v}{d{x}^2}+\rho g=\frac{d{p}_1}{\mathrm{dx}}.$$

(7)

The left side of the equation includes a nonlinear term involving the power $n$ of the velocity gradient and the second derivative of velocity, indicating fluid movement against gravity. In this case, the pressure term connected to gravity will be crucial. Equation (7) is transformed into:

$$a\frac{d^{2}v}{d{x}^2}+\mathrm{nb}{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^{n-1}\frac{d^{2}v}{d{x}^2}-\rho g=0.$$

(8)

The equation incorporates gravity element $\rho g$ to account for the fluid movement against gravity, affecting momentum balance, and incorporating linear and nonlinear velocity gradients to indicate fluid resistance to flow with:

$$v=u_0\mathrm{at}x=0\mathrm{and}{S}_{\mathrm{xy}}=0\mathrm{at}x=\delta .$$

(9)

This outlines the problem’s boundary conditions. The velocity $v$ is set to $u_0$ , the velocity at the origin. The shear stress $S_{\mathrm{xy}}=0\mathrm{at}x=\delta$ , a no-shear condition, meaning that the shear forces disappear at this point, where

$$S_{\mathrm{xy}}=a\frac{\mathrm{dv}}{\mathrm{dx}}+b{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n.$$

(10)

The standard linear shear stress is represented by the first component, $\frac{\mathrm{dv}}{\mathrm{dx}}$ . A nonlinearity, dependent on the velocity gradient increased to the power $n$ , is introduced by the second term, $b{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^n$ . This term describes the non-Newtonian behavior of the fluid, in which the stress exhibits a power-law behavior rather than being proportionate to the rate of strain. Using equation (10) in equation (9) we get:

$$\frac{\mathrm{dv}}{\mathrm{dx}}=0\mathrm{at}x=\delta .$$

(11)

The equation indicates a zero-velocity gradient $\frac{\mathrm{dv}}{\mathrm{dx}}$ at $x=\delta$ , potentially indicating a no-slip or constant velocity at this point, indicating a zero rate of change. This could mean that there is no slippage or that the velocity achieves a constant value and its rate of change zeroes at the boundary $x=\delta$ .

$$\frac{d^{2}v}{d{x}^2}+\frac{\mathrm{nb}}{a}{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^{n-1}\frac{d^{2}v}{d{x}^2}-\frac{\rho g}{a}=0.$$

(12)

The simplified momentum equation combines velocity derivatives with velocity gradient, gravity force per unit area, and a constant, incorporating linear and nonlinear shear stresses, and accounts for the effects on the velocity profile. Substituting $n=3,b=2({\beta }_2+{\beta }_3)$ and $a=\mu$ in equation (12), we have:

$$\frac{d^{2}v}{d{x}^2}+\frac{6({\beta }_2+{\beta }_3)}{\mu }{\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^2\frac{d^{2}v}{d{x}^2}-\frac{\rho g}{\mu }=0.$$

(13)

The equation (12) simplifies to include specific parameters like $n=3,b=2({\beta }_2+{\beta }_3)$ and $a=\mu$ , demonstrating the velocity profile’s dependence on shear stress. A coefficient adjusts the velocity gradient’s influence on the second derivative, while a term $\raisebox{1ex}{$\rho g$}\!\left/ \!\raisebox{-1ex}{$\mu $}\right.$ accounts for gravity’s influence.

$$\begin{array}{c}v=u_0\mathrm{at}x=0\hfill \\ \multicolumn{1}{c}{\frac{\mathrm{dv}}{\mathrm{dx}}=0\mathrm{at}x=\delta }\end{array}\}.$$

(14)

Here, $v^{*}=\frac{v}{u_0},x^{*}=\frac{x}{\delta },b^{*}=6({\beta }_2+{\beta }_3)\frac{u_0^2}{\mu },{\mathrm{\lambda }}^{*}=\frac{\rho g}{\mu u_0}$ are dimensionless parameters.

The problem’s boundary conditions are determined by the velocity gradient at $x=0$ , indicating either flat or steady flow, and zero at $x=\delta$ , indicating inlet velocity. The dimensionless version of (13), considering (14), is expressed without the use of the asterisk symbol as follows:

$$\frac{d^{2}v}{d{x}^2}+\beta {\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^2\frac{d^{2}v}{d{x}^2}-\mathrm{\lambda }=0.$$

(15)

This equation is the dimensionless form of equation (13), obtained by introducing dimensionless parameters. The dimensionless parameter $\beta$ captures the effect of the nonlinear shear stress on the velocity profile, while $g_p$ represents the dimensionless gravity term with:

$$\frac{\mathrm{dv}}{\mathrm{dx}}=0\mathrm{at}x=1\mathrm{and}v=1\mathrm{at}x=0.$$

(16)

The boundary conditions for the dimensionless problem are $x=1$ (normalized boundary), where the velocity gradient is zero, and the dimensionless velocity $v=1$ at $x=0$ . Applying the definitions of fractional calculus, equation (15) can be expressed in fractional form as follows:

$$\frac{d^{2}v\left(x\right)}{d{x}^2}+\beta {\left(\frac{\mathrm{dv}\left(x\right)}{\mathrm{dx}}\right)}^2\frac{d^{2}v\left(x\right)}{d{x}^2}-\mathrm{\lambda }=0.$$

(17)

Equation (17) is a dimensionless form with a non-integer order derivative $D^{\alpha }v\left(x\right)$ , the equation represents an arrangement that permits complicated dynamics in phenomena that display fractional behavior with:

$$v^{\prime }\left(1\right)=0,v\left(0\right)=1.$$

(18)

The formulation’s boundary conditions for fractional calculus define the derivative of velocity as zero at $x=1$ , and $v\left(0\right)=1$ similar to equation (16).

Model for the drainage case

Gravity plays a crucial role in fluid flow, particularly in thin film drainage flow. Gravity accelerates the fluid along the surface, facilitating its downward movement. The gravity-related pressure gradient, a spatial variation of pressure in the fluid, determines the direction and magnitude of fluid flow. In the case of drainage flow, this gradient resists the fluid’s downward motion due to pressure accumulation at the base of the film or the lower parts of the surface.

Equation (15), which originally models thin film flow dynamics with factors like viscosity and pressure gradients, will change when gravity is involved in drainage. The gravitational term and pressure gradient term need to be balanced in the governing equation, accounting for gravity-induced downward motion and the pressure gradient. Thus, gravity facilitates the flow by pulling the fluid downward, but the pressure gradient caused by the fluid’s weight pushes back, affecting the overall flow dynamics.

$$\frac{d^{2}v\left(x\right)}{d{x}^2}+\beta {\left(\frac{\mathrm{dv}\left(x\right)}{\mathrm{dx}}\right)}^2\frac{d^{2}v\left(x\right)}{d{x}^2}+\mathrm{\lambda }=0,$$

(19)

with $v^{\prime }\left(1\right)=0,v\left(0\right)=0$ .

The formulation’s boundary conditions for fractional calculus define the derivative of velocity as zero at $x=1$ , and $v\left(0\right)=0$ .

General framework of Galerkin’s finite element

The following is the procedure used to formulate Galerkin’s finite element procedure, working from the governing differential equation to the matrix formulation:

• I.Convert to standard form: The procedure starts with shifting the governing differential equation to the following form:

  $$f(U,U^{\prime },U^{″})=0$$

where $U$ is the unknown function and the first and second-order spatial derivatives are represented by its derivatives, $U^{\prime }$ and $U^{″}$ , respectively.

• II.Use Quasi-linearization to address the governing equation’s nonlinearity: To transform the governing equation into a linear form, this method requires approximating the nonlinear variables.

  $$f(U_{n+1},U_{n+1}^{\prime },U_{n+1}^{″})\approx \mathrm{linear}\mathrm{terms}\mathrm{in}{U}_n$$

The updated solution from the subsequent iteration is shown above as $U_{n+1}$ . We use the values from the current iteration $U_n$ to linearize the equation about $U_{n+1}$ .

• III.Reduce to a linearized form: To make the equation simpler, the terms about $U_{n+1},U_{n+1}^{\prime }\mathrm{and}{U}_{n+1}^{″}$ are divided and positioned on one side, with the remaining terms on the other:

  $$A_1{U}_{n+1},A_2{U}_{n+1}^{\prime }\mathrm{and}{A}_3{U}_{n+1}^{″}=A,$$

where $A$ , $A_1$ , $A_2$ , and $A_3$ are coefficients that may depend on $U_n,U_n^{\prime }\mathrm{and}{U}_n^{″}$ .

• IV.Create the residual weighted form: The governing equation is linearized and then translated into a weighted residual form. This phase introduces the weighting function $v\left(x\right)$ .

  $$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}v\left(x\right)(A_1{U}_{n+1},A_2{U}_{n+1}^{\prime }\mathrm{and}{A}_3{U}_{n+1}^{″}-A}\end{array})\mathrm{dx}=0.$$

Here, $\mathrm{\Omega }$ is the domain of the problem.

• V.Transform the weighted residual form into a weak form: To do this, the order of the derivatives is lowered by integrating parts (particularly for the second-order derivative $U_{n+1}^{″}$ ):

  $$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}v\left(x\right)U_{n+1}^{″}\mathrm{dx}=-{v\left(x\right)U_{n+1}^{\prime }|}_0^L}\\ \multicolumn{1}{c}{+{\int }_{\mathrm{\Omega }}{v}^{\prime }{U}_{n+1}^{\prime }\mathrm{dx}.}\end{array}$$

The boundary term ${v\left(x\right)U_{n+1}^{\prime }|}_0^L$ appears due to the integration by parts and will be used to impose boundary conditions.

• VI.Split the weak form equation into discrete parts: The weak form equation for the domain $\mathrm{\Omega }$ is stated as a summation over all of its $n$ finite elements. Now, each integral is assessed across every finite element, $e$ .

• VII. Approximate functions using a finite basis: The solution $U_{n+1}^{\prime }\left(x\right)$ and $U_{n+1}\left(x\right)$ are approximated using a set of basis functions ${\psi }_i^{\prime }\left(x\right)$ and ${\psi }_i\left(x\right)$ respectively:

  $$\begin{array}{c}\multicolumn{1}{c}{U_{n+1}^{\prime }\left(x\right)\approx \sum _{i=1}^N{U}_i{\psi }_i^{\prime }\left(x\right),}\\ \multicolumn{1}{c}{U_{n+1}\left(x\right)\approx \sum _{i=1}^N{U}_i{\psi }_i\left(x\right),}\end{array}$$

where $N$ denotes the number of nodes in an element, $U_i$ are the nodal values of the solution, ${\psi }_i\left(x\right)$ are the shape functions associated with the finite element.

• VIII.Apply Galerkin’s finite element approach: It requires the weighting function $v\left(x\right)$ to be chosen as the same set of basis functions used to approximate $U_{n+1}^{\prime }\left(x\right)$ and $U_{n+1}\left(x\right)$ :

  $$v\left(x\right)={\psi }_j\left(x\right).$$

Substituting these approximations into the weak form results in a discretized system of equations.

• IX.Formulate the equation system in matrix notation: A matrix representation of the discretized system of equations is provided. This results in a force vector and a stiffness matrix. The boundary term ${v\left(x\right)U_{n+1}^{\prime }|}_0^L$ enforces appropriate boundary conditions. The final set of equations can be expressed as follows in matrix form:

  $$\mathcal{K}\mathbf{U}=\mathcal{F}$$

This completes Galerkin’s finite element framework resulting in a system of linear equations that can be solved for the unknown nodal values $U$ .

Galerkin’s finite element application to non-Newtonian third-grade Sisko fluid flow

It is necessary to select an appropriate trial or basis function that is applied locally across a typical finite element in the entire domain to use Galerkin’s finite element formulation^44,55 in the finite element environment. We will use $\mathfrak{U}$ to represent this trial function. Compatibility between the elements about displacements must be met in this situation. The trial function is hence $C^0$ -continuous. Two nodes make up each element. At every elemental node, we interpolate the function. Every element node needs to have one unknown parameter for this to work.

$$\mathfrak{U}\left(\chi \right)=b_0+b_1\chi .$$

(20)

The convention employed by Stasa ⁴⁴ and Zienkiewicz ⁴⁵ is to recast the above linear trial function in terms of values of the dependent functions at nodes $i\&j$ , as opposed to expressing the problem in terms of arbitrary constants $b_0$ and $b_1$ .

$$\mathfrak{U}\left(\chi \right)=\left[\begin{array}{cc}{\psi }_1\hfill & {\psi }_2\hfill \end{array}\right]\left[\begin{array}{c}{\mathfrak{U}}_i\hfill \\ \multicolumn{1}{c}{{\mathfrak{U}}_j}\end{array}\right].$$

(21)

The nodal values of the dependent variable $\tilde{\mathfrak{U}}$ are now the trial function constants having the well-known interpolation, or shape, functions ${\psi }_i$ .

Dynamics in lifting case

Considering the lifting case:

$$\frac{d^2\mathfrak{U}\left(\chi \right)}{d{\chi }^2}+\beta {\left(\frac{d\mathfrak{U}}{d\chi }\right)}^2\frac{d^2\mathfrak{U}}{d{\chi }^2}-\lambda =0,$$

(22)

with ${\mathfrak{U}}^{\prime }\left(1\right)=0,\mathfrak{U}\left(0\right)=1$ .

Our systematic methodology will be used to apply quasi-linearization ⁴⁶ to the given nonlinear differential equation. Quasi-linearization is used to replace nonlinear differential equations with linear ones. We employ an iterative process, where:

$$\Phi =\frac{d^2{\mathfrak{U}}_{\mathfrak{n}}}{d{\mathrm{\chi }}^2}+\beta {\left(\frac{d{\mathfrak{U}}_{\mathfrak{n}}}{d\mathrm{\chi }}\right)}^2\frac{d^2{\mathfrak{U}}_{\mathfrak{n}}}{d{\mathrm{\chi }}^2}+\lambda .$$

(23)

The solution’s current approximation is denoted as ${\mathfrak{U}}_n$ . The revised approximation that will be updated is ${\mathfrak{U}}_{n+1}$ . Computing partial derivative of $\Phi$ for ${\mathfrak{U}}_{\mathfrak{n}}^{\prime }$ and ${\mathfrak{U}}_{\mathfrak{n}}^{″}$ :

$$\begin{array}{c}\frac{\partial {\Phi }_n}{\partial {\mathfrak{U}}_{\mathfrak{n}}^{\prime }}=2\beta \frac{d{\mathfrak{U}}_{\mathfrak{n}}}{d\mathrm{\chi }}\frac{d^2{\mathfrak{U}}_{\mathfrak{n}}}{d{\mathrm{\chi }}^2}=2\beta {\mathfrak{U}}_n^{\prime }{\mathfrak{U}}_n^{″},\hfill \\ \multicolumn{1}{c}{\frac{\partial {\Phi }_n}{\partial {\mathfrak{U}}_{\mathfrak{n}}^{″}}=1+\beta {\left(\frac{d{\mathfrak{U}}_{\mathfrak{n}}}{d\mathrm{\chi }}\right)}^2=1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}\end{array}\}.$$

(24)

A common use of the quasi-linearization strategy is to represent the differential equations’ nonlinear part in terms of linear approximations centered on the current estimate, ${\mathfrak{U}}_n$ . This includes: (i) calculating the nonlinear term’s partial derivatives concerning first and second derivatives and (ii) linearizing the terms ${\mathfrak{U}}_{n+1}^{\prime }-{\mathfrak{U}}_n^{\prime }$ and ${\mathfrak{U}}_{n+1}^{″}-{\mathfrak{U}}_n^{″}$ . Partial derivatives of the current approximation ${\mathfrak{U}}_n$ are denoted as ${\mathfrak{U}}_n^{\prime }$ and ${\mathfrak{U}}_n^{″}$ , respectively results in:

$${\left(\frac{\partial \Phi }{\partial {\mathfrak{U}}^{\prime }}\right)}_n({\mathfrak{U}}_{n+1}^{\prime }-{\mathfrak{U}}_n^{\prime })+{\left(\frac{\partial \Phi }{\partial {\mathfrak{U}}^{″}}\right)}_n({\mathfrak{U}}_{n+1}^{″}-{\mathfrak{U}}_n^{″})=0.$$

(25)

Substitute the partial derivatives from (24) into (25) we get:

$$\begin{array}{c}\multicolumn{1}{c}{2\beta \frac{d{\mathfrak{U}}_{\mathfrak{n}}}{d\mathrm{\chi }}\frac{d^2{\mathfrak{U}}_{\mathfrak{n}}}{d{\mathrm{\chi }}^2}({\mathfrak{U}}_{\mathfrak{n}+1}^{\prime }-{\mathfrak{U}}_{\mathfrak{n}}^{\prime })+(1+\beta {\left(\frac{d{\mathfrak{U}}_{\mathfrak{n}}}{d\mathrm{\chi }}\right)}^2)}\\ \multicolumn{1}{c}{({\mathfrak{U}}_{\mathfrak{n}+1}^{″}-{\mathfrak{U}}_{\mathfrak{n}}^{″})=0.}\end{array}$$

(26)

Now substituting ${\mathfrak{U}}_{\mathfrak{n}}^{″}=\frac{\lambda }{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}$ in (26) and rearranging we obtain:

$$\begin{array}{c}\multicolumn{1}{c}{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2){\mathfrak{U}}_{\mathfrak{n}+1}^{″}+2\beta {\mathfrak{U}}_n^{\prime }\frac{\lambda }{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}{\mathfrak{U}}_{\mathfrak{n}+1}^{\prime }}\\ \multicolumn{1}{c}{=\lambda +2\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2\frac{\lambda }{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}.}\end{array}$$

(27)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathfrak{U}}_{n+1}^{″}+2\beta \lambda \frac{{\mathfrak{U}}_n^{\prime }{\mathfrak{U}}_{n+1}^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}=\lambda \frac{1}{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}}\\ \multicolumn{1}{c}{+2\beta \lambda \frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}.}\end{array}$$

(28)

Equation (27) leads to the final linearized form presented in equation (28). For numerical techniques such as finite element formulation, a differential equation can be transformed into its weak form using the weighted residual approach. Here, we’ll obtain the weak formulation of the linearized differential equation by applying the weighted residual technique to it. The differential equation’s residual is denoted by $r\left(\chi \right)$ . Regarding the linearized equation, the residual could be:

$$\begin{array}{c}\multicolumn{1}{c}{r\left(\chi \right)={\mathfrak{U}}_{n+1}^{″}+2\beta \lambda \frac{{\mathfrak{U}}_n^{\prime }{\mathfrak{U}}_{n+1}^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}-\frac{\lambda }{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}}\\ \multicolumn{1}{c}{-2\beta \lambda \frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}.}\end{array}$$

(29)

The residual $r\left(\chi \right)$ defined in equation (29) is multiplied by the Galerkin’s weighting function $w\left(\chi \right)$ , which is often chosen to satisfy the problem’s boundary conditions. Over the domain $\chi ϵ\mathcal{D}=[0,L]$ , integrate the product of the residual and the weight function:

$${\int }_{\mathcal{D}}\mathrm{wr}\left(\chi \right)d\chi =0,$$

(30)

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathcal{D}}w({\mathfrak{U}}_{n+1}^{″}+2\beta \lambda \frac{{\mathfrak{U}}_n^{\prime }{\mathfrak{U}}_{n+1}^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}-\frac{\lambda }{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}}\\ \multicolumn{1}{c}{-2\beta \lambda \frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2})d\chi =0,}\end{array}$$

(31)

Taking the integration by parts of the first term of equation (31), and then rearranging, we will get:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathcal{D}}{w}^{\prime }{\mathfrak{U}}_{n+1}^{\prime }d\chi -2\beta \lambda {\int }_{\mathcal{D}}w{\mathfrak{U}}_{n+1}^{\prime }\frac{{\mathfrak{U}}_n^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}}\\ \multicolumn{1}{c}{d\chi ={w{\mathfrak{U}}_{n+1}^{\prime }|}_0^L-\lambda {\int }_{\mathcal{D}}w\frac{1}{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}d\chi }\\ \multicolumn{1}{c}{-2\beta \lambda {\int }_{\mathcal{D}}w\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .}\end{array}$$

(32)

This is the differential equation’s weak form. It is appropriate for use with finite element methods or other numerical approaches because it only requires the test function $w\left(\chi \right)$ and the first-order derivatives of ${\mathfrak{U}}_{n+1}$ . Equation (32), after the trial functions have been substituted, can be expressed for a specific element ‘ $e_n$ ’ as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{e_n}{w}^{\prime }{\mathfrak{U}}_{n+1}^{\prime }\left(\chi \right)d\chi -2\beta \lambda {\int }_{e_n}w{\mathfrak{U}}_{n+1}^{\prime }\frac{{\mathfrak{U}}_n^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}}\\ \multicolumn{1}{c}{d\chi ={\left[w{\mathfrak{U}}_{n+1}^{\prime }\right]|}_0^L-\lambda {\int }_{e_n}w\frac{1}{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}d\chi }\\ \multicolumn{1}{c}{-2\beta \lambda {\int }_{e_n}w\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .}\end{array}$$

(33)

Approximate ${\mathfrak{U}}_{n+1}^{\prime }\left(\chi \right)\approx \sum _{i=1}^N{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)$ , and $w\left(\chi \right)={\psi }_j\left(\chi \right)$ , where ${\mathfrak{U}}_i$ are the nodal values, ${\psi }_i^{\prime }\left(\chi \right)$ are the derivatives of the shape functions (finite basis functions) for ${\mathfrak{U}}_{n+1}^{\prime }\left(\chi \right)$ , and $N$ is the number of nodes in the element. Replace $w\left(\chi \right)$ and ${\mathfrak{U}}_{n+1}^{\prime }\left(\chi \right)$ with their approximations in the integrals. The weak form becomes:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{e_n}\left(\sum _{i=1}^N{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right){\psi }_j^{\prime }\left(\chi \right)d\chi -2\beta \lambda {\int }_{e_n}\left(\sum _{i=1}^N{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right){\psi }_j}\\ \multicolumn{1}{c}{\left(\chi \right)\frac{{\mathfrak{U}}_n^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi ={\left[{\psi }_j\left(\chi \right){\mathfrak{U}}_{n+1}^{\prime }\right]|}_0^L}\\ \multicolumn{1}{c}{-\lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{1}{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}d\chi }\\ \multicolumn{1}{c}{-2\beta \lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .}\end{array}$$

(34)

Two nodal values for a particular element are given by:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{e_n}\left(\sum _{i=1}^2{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right){\psi }_j^{\prime }\left(\chi \right)d\chi -2\beta \lambda {\int }_{e_n}\left(\sum _{i=1}^2{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right)}\\ \multicolumn{1}{c}{{\psi }_j\left(\chi \right)\frac{{\mathfrak{U}}_n^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi ={\left[{\psi }_j\left(\chi \right){\mathfrak{U}}_{n+1}^{\prime }\right]|}_0^L}\\ \multicolumn{1}{c}{-\lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{1}{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}d\chi -2\beta \lambda {\int }_{e_n}{\psi }_j\left(\chi \right)}\\ \multicolumn{1}{c}{\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .}\end{array}$$

(35)

We can express the above equation in matrix form by collecting terms for each ${\mathfrak{U}}_i$ . Stiffness Matrix ${\mathcal{K}}_{\mathrm{ij}}$ for the first term:

$${\mathcal{K}}_{\mathrm{ij}}^{\left(1\right)}={\int }_{e_n}{\psi }_i^{\prime }\left(\chi \right){\psi }_j^{\prime }\left(\chi \right)d\chi .$$

(36)

Stiffness Matrix ${\mathcal{K}}_{\mathrm{ij}}$ for the second term:

$${\mathcal{K}}_{\mathrm{ij}}^{\left(2\right)}=2\beta \lambda {\int }_{e_n}{\psi }_i^{\prime }\left(\chi \right){\psi }_j\left(\chi \right)\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .$$

(37)

Force Vector ${\mathcal{F}}_i$ for the third term:

$${\mathcal{F}}_i^{\left(1\right)}=\left[\sum _{i=1}^N{\mathfrak{U}}_i{\psi }_j^{\prime }\left(\chi \right){\psi }_j\left(\chi \right)\right]{|}_0^L.$$

(38)

Force Vector ${\mathcal{F}}_i$ for the fourth term:

$${\mathcal{F}}_i^{\left(2\right)}=\lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{1}{1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2}d\chi .$$

(39)

Force Vector ${\mathcal{F}}_i$ for the fifth term:

$${\mathcal{F}}_i^{\left(3\right)}=2\beta \lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .$$

(40)

The global system of equations is assembled by summing contributions from each element over all elements $n$ . This leads to:

$${\mathcal{K}}^{\left(1\right)}\mathfrak{U}-{\mathcal{K}}^{\left(2\right)}\mathfrak{U}={\mathcal{F}}^{\left(1\right)}-{\mathcal{F}}^{\left(2\right)}-{\mathcal{F}}^{\left(3\right)},$$

(41)

where ${\mathcal{K}}^{\left(1\right)}$ and ${\mathcal{K}}^{\left(2\right)}$ are the global stiffness matrices, and ${\mathcal{F}}^{\left(1\right)}$ , ${\mathcal{F}}^{\left(2\right)}$ , and ${\mathcal{F}}^{\left(3\right)}$ are the global force vectors. The final system of equations in matrix notation is:

$$\mathcal{K}\mathfrak{U}=\mathcal{F}.$$

(42)

This matrix equation can then be solved for the unknown nodal values of $\mathfrak{U}$ .

Dynamics in drainage case

When the fluid falls under gravity, the pressure gradient acts in the opposite direction, balancing the system to some extent. This means that the pressure gradient does not assist the flow but opposes it, potentially slowing the flow rate depending on the fluid’s thickness, viscosity, and surface slope. The governing equation (19) for the drainage changes to the following equation (43).

$$\frac{d^2\mathfrak{U}\left(\chi \right)}{d{\chi }^2}+\beta {\left(\frac{d\mathfrak{U}}{d\chi }\right)}^2\frac{d^2\mathfrak{U}}{d{\chi }^2}+\lambda =0,$$

(43)

with ${\mathfrak{U}}^{\prime }\left(1\right)=0,\mathfrak{U}\left(0\right)=0$ .

Following the calculations presented in Subsection “Dynamics in lifting case” we obtain the following relation for a particular element in the drainage case:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{e_n}\left(\sum _{i=1}^2{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right){\psi }_j^{\prime }\left(\chi \right)d\chi +2\beta \lambda {\int }_{e_n}}\\ \multicolumn{1}{c}{\left(\sum _{i=1}^2{\mathfrak{U}}_i{\psi }_i^{\prime }\left(\chi \right)\right){\psi }_j\left(\chi \right)\frac{{\mathfrak{U}}_n^{\prime }}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi }\\ \multicolumn{1}{c}{={\left[{\psi }_j\left(\chi \right){\mathfrak{U}}_{n+1}^{\prime }\right]|}_0^L+\lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{1}{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}d\chi }\\ \multicolumn{1}{c}{+2\beta \lambda {\int }_{e_n}{\psi }_j\left(\chi \right)\frac{{\left({\mathfrak{U}}_n^{\prime }\right)}^2}{{(1+\beta {\left({\mathfrak{U}}_n^{\prime }\right)}^2)}^2}d\chi .}\end{array}$$

(44)

The global system of equations for the drainage case is assembled by summing over all elements similar to the lifting case.

Results and discussions

This article examines the behavior of third-grade non-Newtonian Sisko fluids (NNFs) in thin film under lifting and draining conditions using the finite element method (FEM). The outcomes for various values of the relevant parameters, together with residual errors to prove convergence are shown in Figures 2–13.

Figure 2.

Impact of variations $\beta$ on $\mathfrak{U}\left(\chi \right)$ with $\lambda =1.5$ for lifting case.

Figure 3.

Impact of variation in $\lambda$ on $\mathfrak{U}\left(\chi \right)$ versus $\chi$ with $\beta =1$ for lifting case.

Figure 4.

Impact of simultaneous variations of $\beta$ and $\lambda$ on $\mathfrak{U}\left(\chi \right)$ versus $\chi$ for lifting case.

Figure 5.

Impact of $\beta$ on ${\mathfrak{U}}^{\prime }\left(\chi \right)$ and approximates of ${\mathfrak{U}}^{\prime }\left(0\right)$ for lifting case at $\lambda =1.5$ .

Figure 6.

Impact of $\lambda$ on ${\mathfrak{U}}^{\prime }\left(\chi \right)$ and approximates of ${\mathfrak{U}}^{\prime }\left(0\right)$ for lifting case at $\beta =1$ .

Figure 7.

Mean absolute residuals $\left|r\left(\chi \right)\right|$ on a logarithmic scale for three different lifting cases.

Figure 8.

Impact of variations $\beta$ on $\mathfrak{U}\left(\chi \right)$ with $\lambda =1$ for draining case.

Figure 9.

Impact of variation in $\lambda$ on $\mathfrak{U}\left(\chi \right)$ versus $\chi$ with $\beta =1$ for drainage case.

Figure 10.

Impact of simultaneous variations of $\beta$ and $\lambda$ $\mathfrak{U}\left(\chi \right)$ versus $\chi$ for drainage case.

Figure 11.

Impact of $\beta$ on ${\mathfrak{U}}^{\prime }\left(\chi \right)$ and approximates of ${\mathfrak{U}}^{\prime }\left(0\right)$ for $\lambda =1$ in drainage case.

Figure 12.

Impact of $\lambda$ on ${\mathfrak{U}}^{\prime }\left(\chi \right)$ and approximates of ${\mathfrak{U}}^{\prime }\left(0\right)$ for $\beta =1$ in drainage case.

Figure 13.

Mean absolute residuals $\left|r\left(\chi \right)\right|$ on a logarithmic scale for different drainage cases.

Figure 2 displays how the parameter $\beta$ affects the decay of function $\mathfrak{U}\left(\chi \right)$ . Physically, the parameter β reflects the nonlinearity in the shear-thickening or shear-thinning behavior of the non-Newtonian Sisko fluid. A higher value of β means stronger nonlinear shear effects. Figure 2 shows that the curve flattens out with increasing $\beta$ , suggesting a slower decrease in $\mathfrak{U}\left(\chi \right)$ with increasing $\chi$ . The curve shows a steep decrease in $\mathfrak{U}\left(\chi \right)$ as $\beta$ increases, with a blue dashed curve showing a less steep decline. The steepest curve (black with diamond marks) displays the impact of $\beta =0.1$ and rapidly decreases the velocity profile as $\chi$ grows. Higher value of $\beta$ lowers the rate of reduction, as seen by the less steep decline seen in the blue dashed curve of $\beta =0.3$ when compared to $\beta =0.1$ . The black curve for $\beta =0.5$ shows a more moderate trend of slower decay, while the red solid line for $\beta =0.7$ shows an even slower decrease, maintaining higher values for larger χ. The black dashed curve for $\beta =0.9$ is the flattest, indicating that for high $\beta$ , $\mathfrak{U}\left(\chi \right)$ retains more of its initial value throughout the domain of $\chi$ .

The gravitational parameter λ quantifies the relative strength of gravitational forces compared to viscous forces in the system. Figure 3 shows that the increasing value of λ (stronger gravity) opposes the lifting motion, causing the fluid velocity $\mathfrak{U}\left(\chi \right)$ to decrease faster with distance. A steadier decay is produced by lower $\lambda$ values, whereas a sharper fall is produced by higher $\lambda$ values. The function with larger values of $\lambda$ exhibits an accelerated decay, as seen by the arrow that points in the direction of rising $\lambda$ . All curves exhibit a decreasing trend in $\mathfrak{U}\left(\chi \right)$ , but the rate of decrease varies with $\lambda$ . The black curve with diamond marks displaying the effect of $\lambda =0.2$ exhibits the slowest decline, suggesting that smaller $\lambda$ values lead to a slower decline in $\mathfrak{U}\left(\chi \right)$ . For most of the domain, the blue dashed curve for $\lambda =0.4$ declines at a slightly faster rate than the $\lambda =0.2$ curve. $\mathfrak{U}\left(\chi \right)$ appears to fall more quickly as $\lambda$ grows, according to the black dashed curve with square markers, which shows a further increase in the rate of decrease. In comparison to the lower $\lambda$ values, the red solid curve shows an even faster fall, with $\mathfrak{U}\left(\chi \right)$ declining more quickly. The steepest fall is seen in the black dashed curve with circle markers, indicating that $\mathfrak{U}\left(\chi \right)$ decays at the fastest rate for larger $\lambda$ values.

The effects of several variations of parameters $(\beta ,\lambda )$ on velocity profile vs $\chi$ are shown in Figure 4. As $\chi$ increases from 0 to 1, $\mathfrak{U}\left(\chi \right)$ drops for all combinations of $(\beta ,\lambda )$ . The slowest decline is displayed by the black curve with diamond markers, which suggests that modest values of both $\beta$ and $\lambda$ lead to a steady decline in in velocity profile. In contrast to $\beta =0.2$ , the blue dashed curve shows a somewhat faster drop, indicating that raising both $\beta$ and $\lambda$ speeds the decay. As both parameters climb further, the black dashed-dotted line with square markers shows a faster decline in $\mathfrak{U}\left(\chi \right)$ , indicating a moderate increase in the decay rate. The red solid line exhibits a notably steeper decrease, indicating that decay rates increase with greater values of $\beta$ and $\lambda$ . Out of all the curves, the black dashed curve with circle markers shows the fastest decline, suggesting that the largest values of $\beta$ and $\lambda$ lead to the greatest reduction in $\mathfrak{U}\left(\chi \right)$ .

Graph in Figure 5 demonstrates how the value of shear-thickening/ thinning parameter $\beta$ strongly influences the early rate of change in velocity of the fluid. However, as χ grows, this effect becomes less pronounced, resulting in a more consistent behavior across a range of $\beta$ values close to $\chi =1$ . The derivative of the displacement or velocity profile is denoted by ${\mathfrak{U}}^{\prime }\left(\chi \right)$ , where $\chi$ is a real number between 0 and 1, apparently corresponding to a nondimensional position or a spatial domain. In a lifting flow scenario, this plot probably shows the derivative of the velocity or displacement profile. The parameter $\beta$ may correspond to a physical attribute such as viscosity, stiffness, or nonlinearity that affects the degree of profile change. The system behaves more “smoothly” around the initial point $(\chi =0)$ as $\beta$ grows. In a lifting scenario, this might be viewed as less sensitivity to initial forces or boundary conditions. This could imply that the system is more resilient to initial deformation or velocity changes for larger $\beta$ . According to the initial values, the initial slope ${\mathfrak{U}}^{\prime }\left(0\right)$ gets less negative as $\beta$ grows, suggesting that the system reacts to the beginning situation less strongly. All curves converge towards a similar profile as $\chi$ approaches 1, indicating that ${\beta }^{\prime }s$ influence becomes less significant for big $\chi$ .

The graph in Figure 6 indicates that the gravitational parameter $\lambda$ has a significant impact on the profile’s initial slope; however, as $\chi$ rises, this effect becomes less noticeable, and the profiles converge at all values of $\lambda$ . In this instance, $\lambda$ probably denotes a quantity that describes the intensity or impact of a physical force (like gravity or lifting force). Higher $\lambda$ values result in stronger changes around $\chi =0$ , as they influence the velocity or displacement gradient near the beginning position (small $\chi$ ). An increased initial gradient or force operating on the system is indicated by a greater (more negative) value of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ near $\chi =0$ as $\lambda$ increases. The significance of $\lambda$ , however, decreases as χ rises, and all profiles tend to converge at $\chi =1$ . The initial slope ${\mathfrak{U}}^{\prime }\left(0\right)$ is almost the same for $\lambda =0.8$ and $\lambda =1.0$ , indicating that there might be a limiting behavior as λ rises above a particular threshold.

The mean absolute residuals $\left|r\left(\chi \right)\right|$ for three distinct lifting scenarios are represented as a function of $\chi$ on a logarithmic scale in Figure 7. The horizontal axis displays $\chi$ , a nondimensional spatial variable, and the vertical axis $\log \left(\mathrm{mean}\left|r\left(\chi \right)\right|\right)$ . The graph shows the mean absolute residuals spanning multiple orders of magnitude, from ${10}^{-15}$ to ${10}^{-11}$ , using a logarithmic scale on the $y$ -axis. This makes it possible to compare minor differences in the residuals across the $\chi$ range more effectively. The behavior of lifting case 1 is generally consistent with minor changes over the domain. As a whole a similar tendency, however, with slightly higher variances, particularly around $\chi =0.4$ and $\chi =0.9$ , is seen in lifting instance 2. When compared to the other two examples, lifting case 3 shows more noticeable residual variations, especially in the middle of the domain $\chi =0.3$ to $\chi =0.6$ .

Figure 8 depicts the impact of the resistance generated by nonlinearity shear parameter β on the velocity of the Sisko fluid in drainage case. The slowest increasing curve indicates that the fluid’s characteristics change more gradually for decreasing levels of $\beta$ . In contrast to $\beta =0.3$ , this curve rises more steeply, suggesting that as $\beta$ grows, the rate of change of $\mathfrak{U}\left(\chi \right)$ increases. The curve keeps rising steeper, indicating that the growth in $\mathfrak{U}\left(\chi \right)$ gets more noticeable as $\beta$ gets closer to $1.0$ . This curve shows a quick climb, indicating that when $\beta$ increases, $\beta$ has a greater influence on fluid dynamics, which accelerates the growth in $\mathfrak{U}\left(\chi \right)$ . The curve with the steepest rise out of all of them shows how, at bigger values of $\beta$ , the fluid behavior responds to changes in $\chi$ significantly more quickly.

The fluid’s behavior in this third-grade fluid drainage model is greatly influenced by the gravitational parameter $\lambda$ . Figure 9 shows that the increase in $\lambda$ leads to a sharper, faster increase in $\mathfrak{U}\left(\chi \right)$ , but decreasing values cause a more gradual increase. The graph indicates that the fluid flow gets more noticeable as $\lambda$ rises, and the arrow highlights this tendency. This demonstrates the direct relationship between the parameter $\lambda$ and the velocity or other important property growth over time or space for the fluid. Since the function $\mathfrak{U}\left(\chi \right)$ increases steadily as $\chi$ approaches 1, and starts at $\chi =0$ , it can be inferred that $\mathfrak{U}\left(\chi \right)$ represents a growing quantity, probable velocity, or other flow property in the model. The slowest rate of rise is seen in this curve for $\lambda =0.3$ , suggesting that the fluid flow develops more gradually for smaller values of $\lambda$ . The rise is more significant than for $\lambda =0.3$ , suggesting that the flow accelerates with increasing $\lambda$ . With a greater climb on this curve at $\lambda =0.9$ , it is evident that higher values of $\lambda$ cause the flow characteristics to change more quickly. As $\lambda$ gets closer to $0.9$ , the red curve increases dramatically, indicating a considerably more significant increase in flow rate or velocity. The curve with the steepest rise indicates that the fluid responds quickly to changes in $\lambda$ , as the flow rate peaks at the highest value of $\lambda$ .

Figure 10 shows how the drainage behavior of a third-grade fluid is affected by various combinations of β and $\lambda$ . A slower drop of $\mathfrak{U}\left(\chi \right)$ is caused by smaller values of $\beta$ and $\lambda$ . It suggests that the fluid maintains the flow properties for a longer period. On the other hand, a steeper decrease is produced by greater values of $\beta$ and $\lambda$ , indicating a faster degradation of the fluid’s characteristics over time or space. This downward tendency is highlighted by the arrow, which shows how the rate of decrease grows as the pair $(\beta ,\lambda )$ increases. The slowest rate of decrease is associated with this order pair $(\beta =0.2,\lambda =0.1)$ , indicating that the fluid’s property holds up better over time than the other pairs. The decrease becomes significantly steeper when $(\beta =0.4,\lambda =0.3)$ , $\beta$ , and $\lambda$ rise, indicating a faster reduction in the fluid’s flow characteristics. The reduction in u(χ) steepens further with this pair $(\beta =0.6,\lambda =0.5)$ , suggesting that higher $\beta$ and $\lambda$ values accelerate the fluid’s flow behavior’s decrease. The order pair $(\beta =0.8,\lambda =0.7)$ curve has a more marked decline, indicating that as these values grow, the factors β and λ have a bigger effect on the flow. Here, the sharpest decline is seen for $(\beta =1.0,\lambda =0.9)$ . When $\beta$ and $\lambda$ are at their maximum levels, the flow characteristics are maximally reduced.

For a drainage situation, when $\lambda$  = 1, the curves in Figure 11 show the evolution of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ as a function of $\chi$ for different values of the parameter $\beta$ . The graphs demonstrate that, for all values of $\beta$ , ${\mathfrak{U}}^{\prime }\left(\chi \right)$ consistently lowers with varying beginning points based on $\beta$ . The curve begins at roughly 0.7422 for $\beta$  = 0.3, and it shows a faster drop than the other values of $\beta$ . For $\beta$  = 0.7, the beginning of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ occurs at around 0.6328, after a little slower decrease. The curve for $\beta$  = 1.0 (black solid line with square markers) begins at approximately 0.5857 and has a more gradual slope. With a smoother decrease than the preceding cases, ${\mathfrak{U}}^{\prime }\left(\chi \right)$ begins at roughly 0.5331, for $\beta$  = 1.5 (red solid line with circular marks). The curve starts at around 0.4969 for $\beta$  = 2.0, indicating the slowest rate of fall over the range of $\chi$ .

In conclusion, the curve becomes less steep and the initial value of ${\mathfrak{U}}^{\prime }\left(0\right)$ drops as $\beta$ increases. This suggests that larger values of $\beta$ lead to a slower reduction in ${\mathfrak{U}}^{\prime }\left(\chi \right)$ throughout the range of $\chi$ . The system’s behavior as $\beta$ varies illustrates how various drainage conditions affect the function ${\mathfrak{U}}^{\prime }\left(\chi \right)$ , with steeper declines occurring at lower $\beta$ levels.

Figure 12 shows the evolution of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ as a function of $\chi$ for various values of the parameter $\lambda$ . The curves for different values of $\lambda$ show a trend toward decreasing values as $\chi$ rises from 0 to 1. When $\lambda =0.3$ (black solid line with diamond markers) is used, ${\mathfrak{U}}^{\prime }\left(\chi \right)$ begins at roughly around 0.2608 and progressively decreases as $\chi$ rises. The starting value of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ is around 0.3822 for $\lambda$  = 0.5. The decrease is comparable but slightly steeper. In the case of $\lambda$  = 0.7, the fall in ${\mathfrak{U}}^{\prime }\left(\chi \right)$ begins at roughly 0.4758 and is more gradual than in the other cases. In the case of $\lambda$  = 0.7, the fall in ${\mathfrak{U}}^{\prime }\left(\chi \right)$ begins at roughly 0.4758 and is more gradual than in the other cases. The red solid line with circle markers, $\lambda$  = 0.9, illustrates a smoother drop that commences at approximately 0.5522, which is the highest point on this set of curves. The curve for $\lambda$  = 1.0, starts at 0.5857 and shows the slowest rate of decline. It follows a similar path to that of $\lambda$  = 0.9.

In conclusion, greater values of $\lambda$ lead to a smoother evolution of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ over the range of $\chi$ . This is because, when $\chi$  = 0, the initial value of ${\mathfrak{U}}^{\prime }\left(\chi \right)$ increases and the decay rate become less steep.

Figure 13 displays the logarithmic scale mean absolute residuals $\left|\gamma \left(\chi \right)\right|$ for three distinct drainage scenarios about the parameter $\chi$ . The residuals in all three cases follow a similar pattern, with a sudden increase at first from around ${10}^{-15}$ to ${10}^{-13}$ , and then a continuous increase as $\chi$ approaches 1. In comparison, drainage case 1 has the highest residual values and shows the greatest changes in residual values over the whole $\chi$ range. Similar behavior can be seen in drainage case 2, which has fewer oscillations and slightly smaller residuals. Although drainage case 3 has the smoothest curve and the lowest total residuals, it nevertheless follows a similar pattern. This shows that drainage case 1 has the largest residuals, suggesting that there may be more variability in that situation, whereas drainage instance 3 performs better in terms of decreasing the residuals. Overall, the trend shows that the residuals develop with increasing $\chi$ , showing different patterns for every drainage instance.

The solution reported by Siddiqui et al. ⁷ using HPM for drainage cases is:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathfrak{U}}_A\left(\chi \right)=\lambda (\chi -\frac{1}{2}{\chi }^2)+6\beta {\lambda }^3}\\ \multicolumn{1}{c}{(\frac{1}{12}{\chi }^4-\frac{1}{3}{\chi }^3+\frac{1}{2}{\chi }^2-\frac{1}{3}\chi )+3{\beta }^2{\lambda }^5}\\ \multicolumn{1}{c}{(\frac{1}{3}\chi -\frac{5}{6}{\chi }^2+\frac{10}{9}{\chi }^3-\frac{5}{6}{\chi }^4+\frac{1}{3}{\chi }^5-\frac{1}{18}{\chi }^6).}\end{array}$$

Table 1 presents the comparisons of absolute residual $r_F\left(\chi \right)$ of solution found by the proposed FEM framework with the absolute residual $r_A\left(\chi \right)$ of the above solution reported in Siddiqui et al. ⁷ The residuals of FEM solutions are very small as compared to those of the above solution for several combinations of values of parameters. Some particular instances are presented in Appendix Notations. One can observe that the absolute residuals $r_A\left(\chi \right)$ stay within the range $[{10}^{-06},{10}^{-04}]$ for $\beta =0.1$ and $\lambda =0.1$ , whereas, $r_F\left(\chi \right)$ lies in the range $[{10}^{-17},{10}^{-13}]$ . Similarly, for $\beta =0.7$ and $\lambda =0.2$ the absolute residual $r_A\left(\chi \right)$ cannot get beyond the range $[{10}^{-04},{10}^{-02}]$ , whereas, whereas, $r_F\left(\chi \right)$ achieves the range $[{10}^{-16},{10}^{-13}]$ . In the third instance of $\beta =0.01$ and $\lambda =0.8$ , $r_F\left(\chi \right)$ remains within the range $[{10}^{-16},{10}^{-13}]$ whereas $r_A\left(\chi \right)$ is as low as ${10}^{-04}$ only. These comparisons show that the accuracy level of the FEM solution is much higher than that of the HPM solution reported in Siddiqui et al. ⁷

Table 1.

Comparison of FEM absolute residuals with previous work for drainage case.

$\chi$	$\beta =0.1,\lambda =0.1$		$\beta =0.7,\lambda =0.2$		$\beta =0.01,\lambda =0.8$
	$r_A\left(\chi \right)$ a	$r_F\left(\chi \right)$	$r_A\left(\chi \right)$ a	$r_F\left(\chi \right)$	$r_A\left(\chi \right)$ a	$r_F\left(\chi \right)$
0	5.00E−04	6.59E−17	2.86 E−02	1.48 E−16	2.58 E−02	7.43 E−16
0.1	4.05 E−04	4.26 E−14	2.31 E−02	9.95 E−14	2.08 E−02	2.27 E−13
0.2	3.20 E−04	3.55 E−14	1.82 E−02	1.71 E−13	1.64 E−02	6.82 E−13
0.3	2.45 E−04	5.68 E−14	1.39 E−02	1.23 E−13	1.26 E−02	2.84 E−13
0.4	1.80 E−04	2.13 E−14	1.02 E−02	5.68 E−14	9.24 E−03	4.55 E−13
0.5	1.25 E−04	2.27 E−13	7.05 E−03	1.28 E−13	6.41 E−03	4.55 E−13
0.6	8.00 E−05	9.95 E−14	4.50 E−03	0	4.10 E−03	3.07 E−12
0.7	4.50 E−05	1.71 E−13	2.53 E−03	4.26 E−14	2.31 E−03	2.84 E−13
0.8	2.00 E−05	2.13 E−13	1.12 E−03	2.13 E−13	1.02 E−03	2.84 E−13
0.8	5.00 E−06	9.24 E−14	2.80 E−04	4.41 E−13	2.56 E−04	2.84 E−13
1	0	1.67 E−13	0	2.84 E−14	0	9.66 E−13

Source: ^aSiddiqui et al. ⁷

Significance and novelty of the current study

This work presented a finite element method-based framework for modeling both the lifting and drainage cases of thin-film third-grade non-Newtonian Sisko fluid flows, contributing a more accurate and robust solution method compared to previously published studies. Although several analytical and approximate methods have been used to study thin film flows of non-Newtonian fluids, a comprehensive Galerkin’s FEM analysis for third-grade and Sisko fluids in both lifting and drainage configurations has not been previously reported. The rationale behind opting Galerkin’s FEM is that Galerkin’s FEM has earned a great reputation in computational mechanics by meticulously uniting mathematical sophistication with practical robustness. It delivers high accuracy for nonlinear systems by minimizing error in a variational framework, supporting adaptive refinement, and preserving physics through weak-form conservation—outperforming alternatives like finite differences that lack these foundational advantages. Therefore, the novelty of the current work lies in the fact that, to the best of the authors’ knowledge, this is the first-ever application of Galerkin’s FEM framework to analyze both the lifting and drainage cases of thin-film third-grade non-Newtonian Sisko fluid flows. The proposed finite element method expands its strength over other numerical and semi-analytical methods in following notably important ways.

• Unlike FDM, which require structured grid points and struggles when geometries are irregular, FEM uses flexible meshing due to which it is much more adaptive to complex domains.

• Perturbation-based techniques (like HPM, HAM) have limited generality due to their heavily dependence on well-defined and small algorithmic parameters. In contrast, FEM is capable of handling high nonlinearities without necessitating such assumptions.

• FEM permits the use of higher-order polynomial shape functions that leads to obtaining solutions with higher precision as compared to many low-order traditional schemes.

• The governing equations for third-grade Sisko fluids are highly nonlinear. Complex boundary conditions are accurately handled without relying on small-parameter approximations. The local discretization and quasi-linearization techniques in FEM’s, ensure its stability and reliability in such nonlinear problems.

• Since Sisko fluids contain intricate interactions between gravity and viscous forces, the authors were inspired to develop a framework that goes beyond a conventional solver. The authors wanted quantifiable endorsement of the obtained solutions strongly supported through residual error analysis. The results presented in this work demonstrate that FEM accomplishes residuals as low as ${10}^{-17}$ , whereas traditional HPM stays around ${10}^{-6}$ to ${10}^{-2}$ . This fact reflects FEM’s much better convergence and accuracy of solution.

Conclusions

This work demonstrated the successful application of the finite element method (FEM) for examining thin film flows of non-Newtonian third-grade fluids regarding lifting and draining scenarios. The study depicts intricate fluid dynamics, even with intrinsic nonlinearities, employing a Galerkin-based FEM technique. The findings demonstrate FEM’s ability to efficiently predict velocity profiles and flow properties of the fluid, particularly how gravitational and material parameters influence these variables. In the lifting scenario, gravitational force affects velocity inversely, whereas the material parameter raises it to counterbalance gravitational attraction. In drainage scenarios, gravity enhances the flow rate while the material constant reduces it. Residuals remain very low in both scenarios confirming the accuracy, superior stability and precision of FEM’s over traditional methods. The methodology emphasizes FEM’s strengths in third-grade fluid studies, making it a useful tool in subsequent engineering and industrial applications. As a future direction, we aim to apply the same approach to more complex fractional models.