Novel Numerical Investigations of Some Problems Based on the Darcy–Forchheimer Model and Heat Transfer

Abstract

In this study, we introduced a new type of basis function and subsequently a Chebyshev delta shaped collocation method (CDSC). We then use this method to numerically investigate both the natural convective flow and heat transfer of nanofluids in a vertical rectangular duct on the basis of a Darcy–Brinkman–Forchheimer model and the electroosmosis-modulated Darcy–Forchheimer flow of Casson nanofluid over stretching sheets with Newtonian heating problems. The approximate solution is represented in terms of Chebyshev delta shaped basis functions. Novel error estimates for interpolating polynomials are derived. Computational experiments were carried out to corroborate the theoretical results and to compare the present method with the existing Chebyshev pseudospectral method. To demonstrate our proposed approach, we also compared the numerical solutions with analytic solutions of the Poisson equation. Computer simulations show that the proposed method is computationally cheap, fast, and spectrally accurate and more importantly the obtained approximate solution can easily be used by researchers in this field.

1. Introduction

Heat transfer process have countless industrial applications, including nuclear reactors, electronic devices, chemical reactors, engines, etc. Because of this fact, over the decades, a growing number of scientists are concerned with rapid heat transfer. Researchers are relentlessly studying to enhance the heat transfer rate and thermal conductivity; in tandem, they are also studying to reduce frictional loss, pressure drop, and pumping power for heat transfer fluid (HTF). A new type of HTF is engineered providing improved thermal properties for heat transfer which is called nanofluid. Maxwell in [1] initiated such a study and proposed nanofluid first time. Later, Choi et al. [2] found experimentally that suspensions of a small volume fraction (less than 1%) of nanoparticles improved the thermal conductivity of classical fluid by around two times. Because of efficient heat transfer properties, nanofluids have many applications in cooling technologies and engineering [3,4,5,6,7].

The mixed convection problem for the fully developed unidirectional flow of Newtonian fluid in a vertical channel was initially investigated by Tao [8]. The mixed convection in vertical channels with asymmetric heated walls problem was studied by Hamadah and Wirtz [9], and their results indicated that the buoyancy effect within the bulk flow would cause a flow reversal near the cold wall and enhance the heat transfer performance near the hot wall.

Xu and Pop [10] proposed an analytical solution for the resulting equation from the fully developed mixed convection flow of Al₂O₃–water, Cu–water, and TiO₂–water nanofluids in a vertical channel problem. They indicated that their analytical solution for the opposing flow only exists for a certain region of the Rayleigh number Ra in a physical sense. They also stated that the nanoparticle volume fraction φ plays a key role for improving the heat and mass transfer characteristics of the fluids. More recently, Siddique et al. [11] considered the natural convection nanofluid flow of an upright channel undergoing chemical reaction and heat absorption for five different nanofluids, namely titanium oxide (TiO₂), aluminum oxide (Al₂O₃), copper oxide (CuO), copper (Cu), and silver (Ag) nanofluids, and they found that Ag–water has a higher temperature due to the higher thermal conductivity of Ag particles as compare to the other nanoparticles Cu, TiO₂, Al₂O₃, and CuO, while Al₂O₃–water has a greater velocity than other nanofluids due to the lower density of Al₂O₃. Chen and Chin in [12] considered heat transfer and entropy generation in fully developed mixed convection nanofluid flow in a vertical channel; they used a differential transformation method (DTM) with the assumption of asymmetric heated walls. They showed that the average entropy generation number of the nanofluid is less than that of pure water. Also, Narendran et al. [13,14] reported a 76.44% enhancement in the Nusselt number for fluid passing through a smooth tube equipped with twisted tape inserts. The investigation showed a notable rise in heat recovery efficiency with the elevation of both the mass flow rate and the volume fraction of the nanofluid.

For the mixed convection in a vertical porous channel problem, where flow is modeled by the Brinkman–Forchheimer-extended Darcy equations, Umavathi et al. [15] used the finite-difference method to solve the resulting equations and they found that viscous dissipation enhances the flow reversal in the case of vertical downward flow while it counters the flow in the case of upward flow. For the case of vertical rectangular duct geometry, Umavathi et al. [15] numerically studied the case of the natural convective flow and heat transfer of nanofluids, on the basis of the Darcy–Forchheimer–Brinkman model. They used the finite-difference method and found that the heat transfer is enhanced due to the use of a nanofluid.

Another important problem is the heat transfer in non-Newtonian fluid flow via a stretching sheet. It is well-known that non-linear relationships between the shear stress and rate of shear strain provide us with the existence of non-Newtonian fluid. Non-Newtonian fluids, like Casson fluids, are considered as shear-thinning fluids. In flow velocity at zero they have infinite viscosity, a yielding tension beyond which no fluid flows, and zero viscosity at zero shear rates. There are several materials such as tomato sauce, honey, and blood which are among the well-known liquids whose behavior can be explained by Casson fluid modelling. Casson nanofluid is the fundamental liquid of Casson nanomaterials. There are numerous studies of Casson nanofluid. For example, Ullah et al. [16] examined heat generation/absorption and magnetized Casson nanofluid stream through a nonlinear permeable stretchable cylindrical media having porosity, they found that the magnitude of the friction factor and mass transfer rate are suppressed with increment in Casson parameter, whereas heat transfer rate is found to be intensified. The transport of heat and mass from the surface of a cylinder coated with a catalyst and subject to an impinging flow of a Casson rheological fluid is investigated by Alizadeh et al. [17]. They found that for low values of the Casson parameter and thus strong non-Newtonian behavior, the porous system has a significant tendency towards maintaining a local thermal equilibrium. Regarding Casson fluid, more recently, Farooq et al. [18] noted that the velocity profile diminishes for the rising estimations of magnetic and mixed convection parameters. To learn more about some further recent works in this area, we suggest consulting the following references [19,20,21,22,23,24,25]. Regarding the Darcy–Forchheimer law, Yasir, Muhammad, et al. in [26] considered numerical aggregation for the dissipative flow of a hybrid nanomaterial. In the Darcy–Forchheimer model problem, they found that a rise in the nanoparticle volume fraction boosts the friction drag coefficient and shows a declining trend toward thermal transfer. Jan, Refat Ullah, et al. in [27] considered heat transfer analysis in a hybrid nanocomposite flow in a stretchable convergent/divergent channel in the presence of Darcy–Forchheimer law and Lorentz forces, and they noticed that the hybrid nanocomposite had more dominant features than the nanocomposite for both scenarios of narrowing/expanding channels.

Advances in microfluidic technology have led to the widespread use of electroosmotic flow (OF) [28] in chemical and biomedical fields, such as DNA separation, cell sorting, ion transport, sample separation, and mixing in microfluidic chips, etc. [29]. EOF is a specific phenomenon of electrodynamic flow in micro- or nano-scale channels.

This flow phenomenon provides the advantages of requiring no external mechanical force, consuming low energy, and being simple to operate. Compared to the conventional-scale channels, the flow within micro-scale channels possesses its own special characteristics, such as relative surface roughness, micro-scale effects, slip effects, surface forces and capillary effects, and rapid heat conduction effects, etc. [30]. Recently, Hafez et al. [31] examined the Casson nanofluid’s Darcy–Forchheimer flow across a stretching sheet. Their analyses are based on the viscous and Joule dissipations that the electroosmosis forces (EOF) have on the Casson nanofluid boundary layer. They used a similarity transform to reduce the governing system of nonlinear partial differential equations to a system of nonlinear ordinary differential equations and finally they used the well-known “bvp-4c” solver in MATLAB for the numerical solution of the system of nonlinear ordinary differential equations.

For the studies and applications above, we note here that the authors obtained the exact analytical solution if possible, and they prefer to use the finite-difference technique or DTM and the obtained pointwise solution. But more recently, other methods have also been used in heating problems; for example, Athith, T. S., et al. in [32] used a genetic algorithm and an artificial neural network in a heat exchanger with a partially filled different high-porosity metal foam problem. For other studies, we refer to the works [33,34,35].

Here, we introduced the new delta shaped Chebyshev collocation method (CDSC) and the Chebyshev pseudospectral (CP) method. There are several advantages of using these methods compared with the finite-difference method. In the DSCM or CP approach, the differential equations are solved pointwise in physical space, where the space derivatives are calculated using base functions. But this method is completely different from finite-difference methods because in this method, the approximate solution is represented as a continuous or piecewise continuous function; hence, we know the approximate solution of the problem at every point. Another advantage of the CDSC or CP method is that the approximate analytical solution can be used by other researchers easily. Regarding the DTM method, this method is about Taylor expansion, and we do not have any error bound, but we have provided here an error analysis for the CP method, which can easily be extended to the CDSC method.

The objective of the present paper is to develop a novel numerical method in general, and as an example, we obtain the approximate analytical solution of the following two different problems and present the following contributions:

• We defined a new delta shaped base function and established it with Chebyshev polynomials.
• We used the Chebyshev delta shaped base functions and introduced the new method called the Chebyshev delta shaped collocation method (CDSC) for interpolation problems.
• We then developed a new algorithm and proof of the new theorem for the error bound of the Chebyshev delta shaped collocation method (CDSC) and established the relation with the Chebyshev pseudospectral method.
• We then applied the CDSC and CP methods to find the approximate solution of the governing equation arising from the natural convective flow and heat transfer of nanofluids in a vertical rectangular duct based on a Darcy–Brinkman–Forchheimer model consideration. We showed that the method can easily be adapted for mixed boundary conditions (the Chebyshev tau delta shaped collocation method).
• We also tested the new CDSC and CP methods to find the approximate solution of the governing equation arising from the electroosmotic Darcy–Forchheimer flow of a Casson nanofluid over stretching sheets with a Newtonian heating problem. Again, we showed that the method can easily be adapted for mixed boundary conditions (the Chebyshev tau delta shaped collocation method).

2. Governing Equations

Before introducing the collocation method, let us consider following two different problems arising from fluid flow:

Problem 1. 

Consider a fully developed steady laminar flow of Newtonian nanofluid saturated with a porous matrix in a long vertical rectangular duct as shown in Figure 1. The porous medium is homogeneous and isotropic. The non-Darcian model includes the effect of inertial forces and a consideration of the effect of viscous and Darcy dissipations is incorporated in defining the model. The fluid is assumed to be incompressible, and the effective viscosity is equal to the viscosity of the fluid. The physical properties characterizing the fluid except density and viscosity are assumed to be constant. For density, we have the Boussinesq approximation to be $\rho ={\rho }_0\left(1-\beta \left(T-T_0\right)\right)$. Therefore, continuity equation is automatically satisfied, and the momentum and energy equation reduce to [13]

$${\mu }_{nf}\left(\frac{{\partial }^{2}W}{\partial X^2}+\frac{{\partial }^{2}W}{\partial Y^2}\right)+g{\left(\rho \beta \right)}_{nf}\left(T-T_0\right)-\frac{{\mu }_{nf}}{DA}W-\frac{{\rho }_{nf}{C}_F}{\sqrt{DA}}W\left|W\right|=0,$$

(1)

and

$$\left(\frac{{\partial }^{2}T}{\partial X^2}+\frac{{\partial }^{2}T}{\partial Y^2}\right)+\frac{{\mu }_{nf}}{K_{nf}DA}{W}^2+\frac{C_f{\rho }_{nf}}{K_{nf}\sqrt{DA}}\left|W\right|W^2=0,$$

(2)

where the effective density of the nanofluid ${\rho }_{nf}=\left(1-\varphi \right){\rho }_r+\varphi {\rho }_s$; here, $\varphi$ represents the solid volume fraction of the nanoparticles. Also, the thermal expansion coefficients and effective dynamic viscosity of the nanofluid are given, respectively, by

$${\left(\rho \beta \right)}_{nf}=\left(1-\varphi \right){\left(\rho \beta \right)}_f+\varphi {\left(\rho \beta \right)}_s \mathrm{a}\mathrm{n}\mathrm{d} {\mu }_{nf}=\frac{{\mu }_f}{{\left(1-\varphi \right)}^{2.5}},$$

and $K_{nf}$ is defined as the thermal conductivity of the nanofluid. For a spherical nanoparticle, it can be written as [1]

$$K_{nf}=K_f⌈\frac{K_s+2K_f-2\varphi \left(K_f-K_s\right)}{K_s+2K_f+2\varphi \left(K_f-K_s\right)}⌉,$$

where the above subscripts nf, f, and s, respectively, refer to the thermophysical properties of the nanofluids, base fluid, and nano-solid particles. After using appropriate dimensionless variables [13], Equations (1) and (2) in dimensionless form become

$$\frac{m}{S^2}\frac{{\partial }^{2}w}{\partial x^2}+m\frac{{\partial }^{2}w}{\partial y^2}+Gr{\left(1-\varphi \right)}^{2.5}\left(1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right)\theta -\frac{w}{Da}-I{\left(1-\varphi \right)}^{2.5}\left(1-\varphi +\varphi \frac{{\left(\rho \right)}_s}{{\left(\rho \right)}_f}\right)\left|w\right|w=0,$$

(3)

$$\frac{1}{S^2}\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}+\frac{Br}{Da{\left(1-\varphi \right)}^{2.5}}\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)w^2+IBr\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)K_{nf}\left|w\right|w^2.$$

(4)

Appropriate boundary conditions for Equations (3) and (4) are given by

$$\begin{array}{c}w=0, k\frac{\partial \theta }{\partial y}=Bi\left(\theta +\frac{1}{2}\right) at y=0 \mathrm{f}\mathrm{o}\mathrm{r} 0\le x\le 1, \\ w=0, k\frac{\partial \theta }{\partial y}=Bi\left(\theta -\frac{1}{2}\right) at y=1 \mathrm{f}\mathrm{o}\mathrm{r} 0\le x\le 1, \\ w=0, \frac{\partial \theta }{\partial x}=0 at x=0 \mathrm{a}\mathrm{n}\mathrm{d} 1 \mathrm{f}\mathrm{o}\mathrm{r} 0\le y\le 1. \end{array}$$

(5)

Here, we used a no-slip boundary condition but the method we used can easily be adjusted regarding no-slip boundary conditions.

Problem 2. 

Recently, in 2024, Hafez et al. [31], consider the Electroosmosis-modulated Darcy–Forchheimer flow of Casson nanofluid over stretching sheets in the presence of Newtonian heating problem as shown in Figure 2 (see for example [31]), we solve the same problem by using delta shaped Chebyshev basis function, they first write governing equations as follows

$$\frac{\partial \widehat{u}}{\partial \widehat{x}}+\frac{\partial \widehat{v}}{\partial \widehat{y}}=0,$$

(6)

$$D_e\widehat{u}=\upsilon \left(1+\frac{1}{\alpha }\right)\frac{{\partial }^2\widehat{u}}{\partial {\widehat{y}}^2}+\widehat{\sigma }\left({\widehat{E}}_0{\widehat{B}}_0-{\widehat{B}}_0^2\widehat{u}\right)+{\rho }_e{E}_x-\frac{\upsilon }{DA}\widehat{u}-\frac{F}{\sqrt{DA}}{\widehat{u}}^2,$$

(7)

$$D_e\widehat{T}=\frac{K_{\widehat{T}}}{\rho C_p}\frac{{\partial }^2\widehat{T}}{\partial {\widehat{y}}^2}+{\widehat{\sigma }\left({\widehat{E}}_0{\widehat{B}}_0-{\widehat{B}}_0^2\widehat{u}\right)}^2+\frac{\mu }{\rho C_p}\left(1+\frac{1}{\alpha }\right){\left(\frac{{\partial }^2\widehat{u}}{\partial {\widehat{y}}^2}\right)}^2,$$

(8)

where $D_e=\widehat{u}\frac{\partial }{\partial \widehat{x}}+\widehat{v}\frac{\partial }{\partial \widehat{y}}$, is the first-order differential operator. The following boundary condition is implemented on the flow [27]:

$$\genfrac{}{}{0pt}{}{\widehat{u}={\widehat{u}}_w, \widehat{v}=0,\frac{\partial \widehat{T}}{\partial \widehat{y}}=-h_s\widehat{T}\left(NH\right),}{\widehat{u}\to 0,\widehat{ T}\to {\widehat{T}}_{\infty }.}$$

(9)

After solving the linear Poisson–Boltzmann equation and using the appropriate similarity transform and dimensionless parameters as in [31], we obtain

$$\left(1+\frac{1}{\alpha }\right)f^{‴}-{f^{\prime }}^2+f{f}^{″}+M\left(E_1-f^{\prime }\right)-F_s{f^{\prime }}^2-\frac{1}{Da{R}_{ex}}{f}^{\prime }+U_{hS}{m}_e^2{e}^{-\eta m_e}=0,$$

(10)

$$\theta ”+PrEc\left(1+\frac{1}{\alpha }\right){f^{\prime }}^2+Prf{\theta }^{\prime }+MEc{\left(f^{\prime }-E_1\right)}^2=0,$$

(11)

and appropriate boundary conditions

$$\genfrac{}{}{0pt}{}{f=0,f^{\prime }=1,{\theta }^{\prime }=-\gamma \left(1+\theta \right) at \eta =0,}{f^{\prime }\to 0,\theta \to 0, as \eta \to \infty .}$$

(12)

The following dimensionless parameters are defined above as [31]

$$Pr=\frac{\upsilon }{{\alpha }_p}, M=\frac{\sigma B_0^2}{a\rho }, {Re}_{\widehat{x}}=\frac{{\widehat{u}}_w\left(x\right)x}{\upsilon }, \gamma =h_s\sqrt{\frac{\upsilon }{a}}, E_1=\frac{E_0}{a\widehat{x}{B}_0}, Ec=\frac{{\widehat{u}}_w^2}{C_p{T}_{\infty }}, F_s=\frac{F}{\sqrt{DA}}\widehat{x}.$$

(13)

We note here that in [31], while solving the Poisson–Boltzmann equation, the authors assume Debye–Hückel linearization. This is not necessary. We can easily use our method to obtain the approximate analytical solution of the Poisson–Boltzmann equation without assuming Debye–Hückel linearization. In this way, we can improve the flow problem described in Equation (10).

3. Delta-Shaped Basis Functions

Definition 1. 

Let $m>0$ be an integer and $q\in \left[1, \infty \right]$. The set of all functions in $L_w^q\left[a,b\right]$ such that all distribution derivatives up to the order m are also in $L_w^q\left[a,b\right]$ is called Sobolev space and this set is denoted by $W_w^{p,q}\left[a,b\right]$ as defined by

$$W_w^{p,q}\left[a,b\right]=\left\{V\in L_w^q\left[a,b\right]: {∆}^{\left(\alpha \right)}V\in L_w^q\left[a,b\right], \left|\alpha \right|\le s\right\},$$

(14)

where $w$ is the weight function, when endowed with the norm

$${‖V‖}_{W_w^{s,q}\left[a,b\right]}=\sum _{\left|\mathit{\alpha }\right|\le \mathit{m}}{‖{∆}^{\left(\mathit{\alpha }\right)}V‖}_{L_w^q\left[a,b\right]}.$$

(15)

In the special case $q=2$, then the Sobolev norm becomes $W_w^{s,2}\left[a,b\right]=H_w^s\left[a,b\right]$ and in particular $L_w^2\left[a,b\right]$ is a weighted space used in this study defined by

$$L_w^2\left[a,b\right]=\left\{V:V \mathrm{m}\mathrm{e}\mathrm{s}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e} \mathrm{a}\mathrm{n}\mathrm{d} {‖V‖}_w^2<\infty \right\}.$$

(16)

It is well-known that the Chebyshev expansion of the function $f\left(x\right)\in L_w^2\left[-1,1\right]$ is

$$\tilde{f}\left(x\right)={\displaystyle {\displaystyle \sum _{k=0}^N}}{a}_k{T}_k\left(x\right).$$

(17)

Then, it follows that the Dirac delta function can be written as [36]

$$\delta \left(x-a\right)=\frac{2}{\pi }\frac{1}{\sqrt{1-x^2}}\left(\frac{1}{2}+{\displaystyle \sum _{i=1}^{\infty }}{T}_i\left(a\right)T_i\left(x\right)\right).$$

(18)

Notice that this series diverges at any point in the interval [−1,1]. Since this function is not defined at the end points of the interval, we define delta shaped Chebyshev functions, which are

$$\delta \left(x-a\right)\cong \frac{1}{2}+{\displaystyle \sum _{i=1}^N}{T}_i\left(a\right)T_i\left(x\right).$$

(19)

We see that this approximation oscillates around the axes; hence, we need regularization. There are several kinds of regularization structure to extract the smooth delta shaped function $I\left(x,\zeta \right)$ and regularized delta shaped function as

$$I_{M,\mathsf{\Theta }}\left(x,\zeta \right)=\frac{1}{2}+{\displaystyle \sum _{i=1}^M}{r}_i\left(M,\mathsf{\Theta }\right)T_i\left(\zeta \right)T_i\left(x\right),$$

(20)

where $r_i\left(M,\mathsf{\Theta }\right)$ is the regularization function (this is called the Chebyshev delta shaped basis function), which can be defined as

1.

Lanczos regularization technique

$$r_i\left(M,\mathsf{\Theta }\right)={{\sigma }_i\left(M\right)}^2, {\sigma }_i\left(M\right)=\frac{\sin \left(\nu \left(i,M\right)\right)}{\nu \left(i,M\right)}, \nu \left(i,M\right)=\frac{i\pi }{M+1} ,$$

(21)

where ${\sigma }_i\left(M\right)$ is the sigma factor of Lanczos and is used to overcome the Gibb’s phenomena in the Fourier expansion of non-smooth functions [37].

2.

Riesz regularization scheme [37]

$$r_i\left(M,\mathsf{\Theta }\right)={\left(1-\frac{i^2}{{\left(M+1\right)}^2}\right)}^{\mathsf{\Theta }}$$

(22)

3.

Abel regularization scheme

$$I_{M,\alpha }\left(x,\zeta \right)=\frac{1}{2}+{\displaystyle \sum _{i=1}^M}{r}_i\left(\alpha \right)T_i\left(\zeta \right)T_i\left(x\right), r_i\left(\alpha \right)=e^{-\alpha {\mu }_i^2},$$

(23)

where $r_{M+1}\left(\alpha \right)=\epsilon$ is a small defined value. This summation is well-known as a heat-kernel regularized sum or a generalized Dirichlet series [37]. Regularized delta shaped basis functions which have compact support can be used in interpolation problems as

$$\tilde{f}\left(x\right)={\displaystyle \sum _{j=1}^N}{d_{j}I}_{M,\mathsf{\Theta }}\left(x,{\zeta }_j\right).$$

(24)

We can call this approximation the Chebyshev delta shaped base function. Subsequently, the collocation method can be called the Chebyshev delta shaped collocation method (CDSC), and we will now show the application of this method to some flow and heat convection problems.

4. Delta Shaped Chebyshev Collocation Method

We now expand the unknown functions (3) and (4) in terms of the delta shaped base function as $\left(s=4,8,12 \mathrm{u}\mathrm{s}\mathrm{e}\mathrm{d}\right)$

$$I_{L}w\left(x,y\right)={\displaystyle \sum _{i=0}^L}{\displaystyle \sum _{j=0}^L}{a}_{ij}\left\{\left[\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(x\right)T_n\left(x\left[i\right]\right)\right)\right]\left[\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[j\right]\right)\right)\right]\right\},$$

(25)

and

$$I_L\theta \left(x,y\right)={\displaystyle \sum _{i=0}^L}{\displaystyle \sum _{j=0}^L}{b}_{ij}\left\{\left[\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(x\right)T_n\left(x\left[i\right]\right)\right)\right]\left[\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[j\right]\right)\right)\right]\right\},$$

(26)

$$A={\left({\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n^{″}\left(x_j\right)T_n\left(x\left[i\right]\right)\right)}_{j=1}^{L-1},$$

(27)

$$B={\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(x_j\right)T_n\left(x\left[i\right]\right)\right)}_{j=1}^{L-1},$$

(28)

and substitute into (3) and (4); thus, we obtain

$$\left(\frac{m}{S^2}A\otimes B+mB\otimes A-\frac{1}{Da}B\otimes B-I{\left(1-\varphi \right)}^{2.5}\left(1-\varphi +\varphi \frac{{\left(\rho \right)}_s}{{\left(\rho \right)}_f}\right)\left|\left(B⨂B\right){\Xi }^n\right|\right){\Xi }^{n+1}=-Gr{\left(1-\varphi \right)}^{2.5}\left(1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right)\left(B\otimes B\right){\mathsf{\Theta }}^n,$$

(29)

and

$$\left(\frac{1}{S^2}A\otimes B+B\otimes A\right){\mathsf{\Theta }}^{n+1}=-\frac{{\mu }_{nf}}{K_{nf}DA}{\left(\left(B\otimes B\right){\Xi }^{n+1}\right)}^2-\frac{C_f{\rho }_{nf}}{K_{nf}\sqrt{DA}}\left|\left(B⨂B\right){\Xi }^{n+1}\right|{\left(\left(B\otimes B\right){\Xi }^{n+1}\right)}^2 ,$$

(30)

where ${\Xi }^{n+1}$ is the $\left(L-1\right)(L-1)$-vector:

$${\Xi }^{n+1}={\left[a_{11},\dots ,a_{1L-1}\left| a_{21},\dots ,a_{2L-1}\right|\dots |a_{\left(L-1\right)1},\dots ,a_{\left(L-1\right)\left(L-1\right)}\right]}^T,$$

(31)

and ${\mathsf{\Theta }}^{n+1}$ is the $\left(L-1\right)(L-1)$-vector:

$${\mathsf{\Theta }}^{n+1}={\left[b_{11},\dots ,b_{1L-1}\left| b_{21},\dots ,b_{2L-1}\right|\dots |b_{\left(L-1\right)1},\dots ,b_{\left(L-1\right)\left(L-1\right)}\right]}^T.$$

(32)

Boundary conditions can be written as $\left(k=0\right)$

$$\left(B\left(\pm 1\right)\otimes B\right){\Xi }^{n+1}=0,\left(B\otimes B\left(\pm 1\right)\right){\Xi }^{n+1}=0,$$

(33)

$$\left(A1\left(\pm 1\right)\otimes B\right){\mathsf{\Theta }}^{n+1}=0, \left(B\otimes B\left(\pm 1\right)\right){\mathsf{\Theta }}^{n+1}=\pm \frac{1}{2},$$

(34)

where $B(x)={\left(\frac{1}{2}+{\displaystyle {\sum }_{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(x\right)T_n\left(x\left[i\right]\right)\right)}_{i=1}^{L-1}$ and $A1\left(x\right)={\left({\displaystyle {\sum }_{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n^{\prime }\left(x\right)T_n\left(x\left[i\right]\right)\right)}_{i=1}^{L-1}$. In this way, we obtain a system of linear equations. We used QR decomposition to find the inverse which is an approximate solution, and this approximate solution is repeated iteratively until the differences between two successive iterations are less than ${10}^{-5}$.

To apply the method for Problem 2, we need the transform of the independent variables from $\left[0,S\right]\to \left[-1,1\right]$. After, we introduce a new dependent variable as $f^{\prime }=g$. Then, we obtain from (10)

$$\frac{8}{S^3}\left(1+\frac{1}{\alpha }\right)g^{″}-\frac{4}{S^2}\left(1+F_s\right)g^2+\frac{4}{S^2}{g}^{\prime }{\displaystyle \underset{-1}{\overset{y}{\int }}}g\left(z\right)dz+M\left(E_1-\frac{2}{S}g\right)-\frac{2}{SDa{R}_{ex}}g+U_{hS}{m}_e^2{e}^{-\left(\frac{yS+S}{2}\right)m_e}=0,$$

(35)

$$g\left(-1\right)=\frac{S}{2}, g\left(1\right)=0,$$

(36)

and

$$\frac{4}{S^2}\theta ”+\frac{4}{S^2}PrEc\left(1+\frac{1}{\alpha }\right)g^2+\frac{2}{S}Pr{\displaystyle \underset{-1}{\overset{y}{\int }}}g\left(z\right)dz{\theta }^{\prime }+MEc{\left(\frac{2}{S}g-E_1\right)}^2=0,$$

(37)

$${\theta }^{\prime }=-\frac{S\gamma }{2}\left(1+\theta \right) at y=-1, \theta \left(1\right)=0.$$

(38)

Now, we assume the solution in the form

$$I_{L}g={\displaystyle \sum _{i=0}^L}{a}_i\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[i\right]\right)\right),$$

(39)

$$I_L\mathsf{\Theta }={\displaystyle \sum _{i=0}^L}{b}_i\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[i\right]\right)\right),$$

(40)

$$\frac{8}{S^3}\left(1+\frac{1}{\alpha }\right)A{\zeta }^{n+1}-\frac{4}{S^2}\left(1+F_s\right)\left(B{\zeta }^n\right)B{\zeta }^{n+1}+\frac{4}{S^2}C{\zeta }^{n+1}{\displaystyle \underset{-1}{\overset{y}{\int }}}B\left(z\right){\zeta }^{n}dz-\frac{2}{S}\left(M+\frac{1}{Da{R}_{ex}}\right) B{\zeta }^{n+1}=-M{E}_1-U_{hS}{m}_e^2{e}^{-\left(\frac{yS+S}{2}\right)m_e},$$

(41)

$$B\left(-1\right){\zeta }^{n+1}=\frac{S}{2}, B\left(1\right){\zeta }^{n+1}=0.$$

(42)

$$\frac{4}{S^2}A{\xi }^{n+1}+\frac{4}{S^2}PrEc\left(1+\frac{1}{\alpha }\right){\left(B{\zeta }^{n+1}\right)}^2+\frac{2Pr}{S}C{\xi }^{n+1}{\displaystyle \underset{-1}{\overset{y}{\int }}}B{\zeta }^{n+1}dz+\frac{4}{S^2}MEc{\left(B{\zeta }^{n+1}\right)}^2-\frac{MEc{E}_1}{S}B{\zeta }^{n+1}=-MEc{E}_1^2,$$

(43)

with

$$A1{\left(-1\right)\xi }^{n+1}+\frac{S\gamma }{2}B\left(-1\right){\xi }^{n+1}=-\frac{S\gamma }{2}.$$

(44)

These algebraic equations were solved iteratively and the result documented in the last section. We also used the Chebyshev pseudospectral method for the approximate solution of Equations (3), (4), (10), and (11).

Problem 3. 

In this problem, we consider the following Poisson equation:

$$-\mathsf{\Delta }u=1, \mathsf{\Omega }=\left(-1,1\right)x\left(-1,1\right), {\left.u\right|}_{\partial \mathsf{\Omega }}=0$$

(45)

with an exact solution,

$$u\left(x,y\right)=-\frac{64}{{\pi }^4}{\displaystyle \sum _{\begin{array}{c}n,m=1\\ n,m odd\end{array}}^{\infty }}{\left(-1\right)}^{\frac{n+m}{2}}\frac{\cos \left(\frac{n\pi x}{2}\right)\cos \left(\frac{m\pi y}{2}\right)}{nm\left(n^2+m^2\right)},$$

(46)

This is the used benchmark problem in a famous paper from Jie Shen in [38] and D. B. Hoidvogel and T.A. Zang [39]. We solved the same problem with the Chebyshev delta shaped collocation method. In

Table 1. Maximum pointwise error $u-u_n$.

N	LGM	CTN	CCM	CDSC	CTDSC	FD2	FD4
16	1.42 × 10−6	3.52 × 10−5	7.47 × 10−7	8.5 × 10−6	3.5301 × 10−5	7.17 × 10−6	9.02 × 10−4

, we list the maximum pointwise error of LGM (Legendre–Galerkin method), CTM (Chebyshev-tau), CCM (Chebyshev collocation method introduced by Shen), CDSC, CTDSC (Chebyshev tau delta shaped collocation method), and the second- and fourth-order finite-difference (FD2 and FD4) methods.

From the table, we see that the method used by Shen [38] gives a better convergence rate, but this method is restricted on zero Dirichlet boundary conditions. The method we used is very flexible with regards to boundary conditions.

Conjecture 1. 

Is there any relationship between the delta shaped Chebyshev bases collocation method and the Chebyshev pseudospectral method? The answer is yes, this relation can be described by matrix transformation. The proof can be seen in the next section.

Conjecture 2. 

If we use same collocation points in both the delta shaped Chebyshev collocation and Chebyshev pseudospectral method, do we obtain same or approximately same solution from both method? The answer depends on the bound of the differences between the solutions (error bound); to show this, we use ten Chebyshev Gauss collocation points, and the difference between the solution and derivative of the solution can be seen in Figure 3 and Figure 4. We see that the error bound is less than $1.2\times {10}^{-3}$. But for $B{\zeta }^n=0$, Equation (43) becomes linear such that an exact analytical solution exists. In this case, given the aforementioned prediction of the two methods compared with the exact analytical solution, we found that our method approached the exact solution faster that Chebyshev pseudospectral method. This shows that the method we describe here is competitive and easy to use.

5. Error Analysis for New Base Functions

In this section, we find the error bound for the numerical solution. To do this, we first look at the expansion

$$\begin{gathered} I_{L}g={\displaystyle \sum _{i=0}^L}{a}_i\left(\frac{T_0\left(y\right)}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[i\right]\right)\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=0.5\left(a_0{T}_0\left(y_0\right)+a_1{T}_0\left(y_2\right)+\cdots +a_L{T}_0\left(y_L\right)\right)T_0\left(y\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+0.99096\left(a_0{T}_1\left(y_0\right)+a_1{T}_1\left(y_2\right)+\cdots +a_L{T}_1\left(y_L\right)\right)T_1\left(y\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+0.357489\left(a_0{T}_{10}\left(y_0\right)+a_2{T}_{10}\left(y_2\right)+\cdots +a_L{T}_{10}\left(y_L\right)\right)T_{10}\left(y\right)+\cdots \end{gathered}$$

(47)

If we remember, the expansion of the function in terms of Chebyshev polynomials is

$$I_{L}g={\displaystyle \sum _{i=0}^L}{b}_{i }{T}_i\left(y\right),=b_0{T}_0\left(y\right)+b_1{T}_1\left(y\right)+\cdots .$$

(48)

Comparing (45) with (46), the two expansions are equal if the coefficient of $T_j\left(y\right),j=0,1\cdots L$, i.e.,

$$\begin{array}{c}\begin{array}{c}0.5\left(a_0{T}_0\left(y_0\right)+a_1{T}_0\left(y_2\right)+\cdots +a_L{T}_0\left(y_L\right)\right)=b_0,\\ 0.99096\left(a_0{T}_1\left(y_0\right)+a_1{T}_1\left(y_2\right)+\cdots +a_L{T}_1\left(y_L\right)\right)=b_{1,}\\ --\end{array}\\ A_L\left(a_0{T}_K\left(y_0\right)+a_1{T}_K\left(y_2\right)+\cdots +a_L{T}_K\left(y_L\right)\right)=b_L.\end{array}$$

(49)

Regarding the solution of (49) which can be written as $T\mathit{a}=\mathit{b}$, the condition number of the coefficient matrix $T$ depends on the selection of the collocation points, where we selected the Chebyshev point as a collocation point, then we found that condition number relatively large. This may be one of the reasons why the two approaches do not give the same approximate solution. Since the determinant of the coefficient matrix $T$ is not zero, then there exists a unique solution of (49); then, if we provide an error analysis for the pseudospectral Chebyshev method, the same error bound is valid for the delta shaped method.

Theorem 1. 

If there is unique solution of the system of equations in (47) and$g\left(y\right)\in H_w^s\left[-1,1\right]$ with $m>0$, then there exists a constant $C$, such that the following error bound holds:

$${‖g-I_{L}g‖}_{L_w^2\left[-1,1\right]}\le C L^{-s}{‖g‖}_{H_w^s\left[-1,1\right]}.$$

(50)

Proof. 

We now approximate the solution obtained from truncated delta shaped Chebyshev base functions as

$$P_{L}g={\displaystyle \sum _{i=0}^L}{a}_i\left(\frac{1}{2}+{\displaystyle \sum _{n=1}^K}{\left(1-\frac{n^2}{{\left(K+1\right)}^2}\right)}^s{T}_n\left(y\right)T_n\left(y\left[i\right]\right)\right),$$

(51)

for $s=4$, we assumed the solution system in (49) exists and is unique. Then, (51) can be written as

$$P_{L}g={\displaystyle \sum _{i=0}^L}{b}_{i }{T}_i\left(y\right),$$

(52)

and we then obtain the expansion coefficient as

$$b_{i }=\frac{1}{{\rho }_i}{\displaystyle \underset{-1}{\overset{1}{\int }}}{P}_{L}g{T}_i\left(y\right)w\left(y\right)dy,{\rho }_i\mathrm{i}\mathrm{s} \mathrm{a} \mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}.$$

(53)

Since from [40], we have

$$n^2{T}_n\left(y\right)w\left(y\right)={\mathcal{L}}_n{T}_n\left(y\right)w\left(y\right)=\frac{\partial }{\partial y}\left(\left(-\sqrt{1-y^2}\right)\frac{\partial }{\partial y} T_n\left(y\right) \right),$$

(54)

using (51) in (50), we obtain

$$b_{i }=\frac{1}{{i^2\rho }_i}{\displaystyle \underset{-1}{\overset{1}{\int }}}{P}_{L}g\frac{\partial }{\partial y}\left(\left(-\sqrt{1-y^2}\right)\frac{\partial }{\partial y} T_i\left(y\right) \right)dy.$$

(55)

Applying integration by parts twice, we have

$$b_{i }=\frac{1}{{i^2\rho }_i}{\displaystyle \underset{-1}{\overset{1}{\int }}}\mathcal{L}\left(P_{L}g\right)T_i\left(y\right)w\left(y\right)dy$$

(56)

which can be applied $m$ times as to obtain

$$b_{i }=\frac{1}{{i^{2m}\rho }_i}{\displaystyle \underset{-1}{\overset{1}{\int }}}{\mathcal{L}}^m\left(P_{L}g\right)T_i\left(y\right)w\left(y\right)dy.$$

(57)

Then, the Cauchy–Schwarz inequality yields the upper error bound for coefficient as

$${\left|b_{i }\right|}^2\le \frac{C}{i^{4m}}{‖{\mathcal{L}}^m\left(P_{L}g\right)‖}_{L_w^1\left[-1,1\right]}^2\le \frac{C}{i^{4m}}{‖P_{L}g‖}_{H_w^{2m}\left[-1,1\right]}^2.$$

(58)

Now, we consider the discrete approximation

$$I_{L}g={\displaystyle \sum _{i=0}^L}{\mathcal{F}}_i\left(y\right)g\left(y_i\right)={\displaystyle \sum _{i=0}^L}{T}_i\left(y\right){\pi }_i,$$

(59)

where ${\pi }_i$ is the discrete expansion coefficient and is defined by

$${\pi }_i=\frac{1}{{\rho }_i}{\displaystyle \sum _{i=0}^L}{T}_i\left(y_i\right)w_{i}g\left(y_i\right).$$

(60)

Under the assumption of sufficient smoothness functions, the aliasing error is

$${\pi }_i=b_i+\frac{1}{{\rho }_i}{\displaystyle \sum _{j>L}^{\infty }}{\left[T_i\left(x\right),T_j\left(x\right)\right]}_w{b}_j,$$

(61)

Here, ${\left[.,.\right]}_w$ is the weighted discrete inner product. Since

$$g\left(y\right)={\displaystyle \sum _{i=0}^{\infty }}{b}_i{T}_i\left(y\right),$$

(62)

using (56)–(59), it is easy to see that

$${‖g-I_{L}g‖}_{L_w^2\left[-1,1\right]}^2={‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}^2+{‖A_{L}g‖}_{L_w^2\left[-1,1\right]}^2,$$

(63)

where

$$A_{L}g={\displaystyle \sum _{i=0}^L}\frac{1}{{\rho }_i}{\displaystyle \sum _{j>L}^{\infty }}\left[T_i\left(y\right),T_j\left(y\right)\right]T_i\left(y\right)b_i={\displaystyle \sum _{j>L}^{\infty }}{\displaystyle \sum _{i=0}^L}\frac{1}{{\rho }_i}\left[T_i\left(y\right),T_j\left(y\right)\right]T_i\left(y\right)b_i.$$

(64)

since $T_i\left(y\right),T_j\left(y\right)$ are orthogonal and the range of summation is different; hence, we obtain $A_{L}g=0$. Therefore,

$${‖g-I_{L}g‖}_{L_w^2\left[-1,1\right]}^2={‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}^2.$$

(65)

On the other hand,

$$g-P_{L}g={\displaystyle \sum _{i>L}^{\infty }}{b_{i}T}_i\left(y\right),$$

(66)

then, we have

$${‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}^2={\displaystyle \sum _{i>L}^{\infty }}{\displaystyle \underset{-1}{\overset{1}{\int }}}{\left|{b_{i}T}_i\left(y\right)\right|}^{2}w\left(y\right)dy={\displaystyle \sum _{i>L}^{\infty }}{\left[T_i\left(y\right),T_i\left(y\right)\right]}_w{\left|b_i\right|}^2={\displaystyle \sum _{i>L}^{\infty }}{\rho }_i{\left|b_i\right|}^2,$$

(67)

using (58), we obtain

$${‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}^2\le C{‖g‖}_{H_w^{2m}\left[-1,1\right]}^2{\displaystyle \sum _{i>L}^{\infty }}{\rho }_i{i}^{-4m}\le C{L}^{-4m}{‖g‖}_{H_w^{2m}\left[-1,1\right]}^2,$$

(68)

or

$${‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}\le C{L}^{-2m}{‖g‖}_{H_w^{2m}\left[-1,1\right]}.$$

(69)

If we say $s=2m$, we obtain

$${‖g-P_{L}g‖}_{L_w^2\left[-1,1\right]}\le C{L}^{-s}{‖g‖}_{H_w^s\left[-1,1\right]} \mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o} {‖g-I_{L}g‖}_{L_w^2\left[-1,1\right]}\le C{L}^{-s}{‖g‖}_{H_w^s\left[-1,1\right]}.$$

(70)

This completes the proof. □

Further, we discuss the bounds of the given electroosmosis-modulated Darcy–Forchheimer flow of a Casson nanofluid over stretching sheets with Newtonian heating problems. We consider the error function $E\left(y\right)=\left(g\left(y\right)-I_{L}g\left(y\right)\right)$, where $g\left(y\right)$ is the exact solution and $I_{L}g\left(y\right)$ is the approximate solution of a given problem obtained from the pseudospectral Chebyshev method. Then, from (35), we have

$$R\left(y\right)=\frac{8}{S^3}\left(1+\frac{1}{\alpha }\right)\left(g^{″}-I_L{g}^{″}\right)-\frac{4}{S^2}\left(1+F_s\right)\left(g^2-I_L{g}^2\right)+\frac{4}{S^2}{g}^{\prime }{\displaystyle \underset{-1}{\overset{y}{\int }}}g\left(z\right)dz-\frac{4}{S^2}{I}_L\left(g^{\prime }\right){\displaystyle \underset{-1}{\overset{y}{\int }}}{I}_L\left(g\left(z\right)\right)dz-\left(\frac{2}{S}+\frac{2}{SDa{R}_{ex}}\right)\left({g-I}_{L}g\right)=0.$$

(71)

since $\left(g^2-I_L{g}^2\right)=\left(g-I_{L}g\right)\left(g+I_{L}g\right)$ and $\frac{4}{S^2}{g}^{\prime }{\displaystyle {\int }_{-1}^y}g\left(z\right)dz-\frac{4}{S^2}{I}_L\left(g^{\prime }\right){\displaystyle {\int }_{-1}^y}{I}_L\left(g\left(z\right)\right)dz=\frac{4}{S^2}\left(g^{\prime }-I_L\left(g^{\prime }\right)\right){\displaystyle {\int }_{-1}^y}g\left(z\right)dz+\frac{4}{S^2}{I}_L\left(g^{\prime }\right){\displaystyle {\int }_{-1}^y}\left({g\left(z\right)-I}_L\left(g\left(z\right)\right)\right)dz$. Taking the $L_w^2\left[-1,1\right]$ norm of both sides of (68), we obtain

$${‖R\left(y\right)‖}_{L_w^2\left[-1,1\right]}^2\le C\left[{‖g^{″}-I_L{g}^{″}‖}_{L_w^2\left[-1,1\right]}^2+{‖g-I_{L}g‖}_{L_w^2\left[-1,1\right]}^2+{‖g^{\prime }-I_L{g}^{\prime }‖}_{L_w^2\left[-1,1\right]}^{2}max{\left({\displaystyle \underset{-1}{\overset{y}{\int }}}g\left(z\right)dz\right)}^2+{‖I_L{g}^{\prime }‖}_{L_w^2\left[-1,1\right]}^{2}max{\left({\displaystyle \underset{-1}{\overset{y}{\int }}}\left({g\left(z\right)-I}_L\left(g\left(z\right)\right)\right)dz\right)}^2\right].$$

(72)

Since $\left|{g\left(z\right)-I}_L\left(g\left(z\right)\right)\right|={\displaystyle {\sum }_{i>L}^{\mathrm{\infty }}}\left|{b_{i}T}_i\left(z\right)\right|$ and $\left|b_i\right|\le \frac{C}{i^{2m}}{‖g‖}_{H_w^{2m}\left[-1,1\right]}$ and using (68), we obtain

$${‖R\left(y\right)‖}_{L_w^2\left[-1,1\right]}^2\le C\left[L^{-4m}{‖g‖}_{H_w^{2m+2}\left[-1,1\right]}^2+L^{-4m}{‖g‖}_{H_w^{2m+1}\left[-1,1\right]}^2+L^{-4m}{‖g‖}_{H_w^{2m}\left[-1,1\right]}^2\right],$$

(73)

which show the method, we used here converges spectrally.

6. Results and Discussion

The representative approximate solution obtained from the new method is displayed in this section graphically for two different values of the Darcy number $Da$, while the other parameters are fixed. However, the method we developed can easily be adopted for other cases. Let us begin with the Darcy number $Gr=10 \mathrm{a}\mathrm{n}\mathrm{d} Da=0.1$ for the first problem. Figure 5a,b represent the upward and downward velocity contours, and these upward and downward velocity distributions are not symmetrical. It is clear from the momentum equations that an increase in the value of the Darcy number causes a decrease in the value of the velocity. Next, we check the temperature distribution corresponding to velocity distribution for $Gr=10\mathrm{a}\mathrm{n}\mathrm{d} Da=0.1$. We see nonsymmetrical contours with respect to the x and y variables, and the corresponding three-dimensional temperature field is given by Figure 5c,d. From the graph in Figure 5c,d, one cannot distinguish the changing temperature distribution with respect to y for fixed x; to see the variation in the velocity and temperature profiles, we provide Figure 6 and Figure 7, which show that both solutions are not symmetrical with respect to y and vary with y for fixed x. In these solutions, we used 28 base functions, means that 1568 unknown coefficients are calculated iteratively, and the stopping criteria here are $\left|w_{i+1}\left(x,y\right)-w_i\left(x,y\right)\right|\le {10}^{-5}$ and $\left|{\theta }_{i+1}\left(x,y\right)-{\theta }_i\left(x,y\right)\right|\le {10}^{-5}$. When we increase the value of the Darcy number, we find that this causes a decrease in the value of velocity and this effects a decrease in the variations in temperature; these are shown in Figure 8d. Next, we used the delta shaped Chebyshev base function for the solution for the electroosmosis-modulated Darcy–Forchheimer flow of a Casson nanofluid over stretching sheets with Newtonian heating problems. For this problem, we assumed that the flow changes happen to be on a finite interval outside of a finite interval velocity of almost zero. The most important parameters mathematically in the equation affect the electroosmotic $\left(U_{hS}\right)$ and electromagnetohydrodynamic $\left(M\right)$ forces on the velocity $f^{\prime }\left(y\right)$, because this directly affects the decay of the solution toward zero. Let us consider both $U_{hS} \mathrm{a}\mathrm{n}\mathrm{d} M=0$; then, we check the effect of the Darcy number on the velocity field, which can be seen in Figure 9. We see that $f^{\prime }\left(y\right)$ quickly decays to zero; but for non-zero values of $U_{hS} \mathrm{a}\mathrm{n}\mathrm{d} M=0$, this is not the same case, as shown in Figure 10. Finally, Figure 11 shows the temperature variation for different values of the electroosmotic $\left(U_{hS}\right)$ and electromagnetohydrodynamic $\left(M\right)$ forces. We also compared our result with the prediction of the Chebyshev pseudospectral method. We see a slightly different result for a small number of base functions; when we increase the base function, we see that the difference between the solutions is negligible, less then ${10}^{-9}$.

7. Conclusions

In this research, we first describe the Chebyshev delta shaped collocation method and use this method to find the approximate analytical solution to two benchmark heat transfer problems which arise from fluid flow problems.

• There are several papers concerning the delta shaped sine collocation method (see for example [41,42,43]), but it is well-known that the Chebyshev pseudospectral method gives a much more accurate solution than the Fourier pseudospectral method for problems with non-periodic boundary conditions. This is why we introduced new Chebyshev delta shaped basis functions.
• By using the Chebyshev delta shaped collocation method, we have solved two benchmark heat transfer problems, and our approximate solution is also original. Since all the other works have used the finite-difference method, their solution is only valid at discrete points.
• Since we use delta shaped base functions with compact support, in this regard, this method can be regarded as between a spectral element method and a Galerkin spectral method.
• A more general error bound is necessary without using the Chebyshev pseudospectral method; such a study under consideration.