# Wavelet Operational Matrices and Lagrange Interpolation Differential Quadrature-Based Numerical Algorithms for Simulation of Nanofluid in Porous Channel

Authors: Aisha M. Alqahtani; Ram Jiwari
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This work analyses the features of nanofluid flow and thermal transmission (NFTT) in a rectangular channel which is asymmetric by developing two numerical algorithms based on scale-2 Haar wavelets (S2HWs), Lagrange’s interpolation differential quadrature technique (LIDQT), and quasilinearization process (QP). In the simulation procedure, first of all, using similarity transformation (ST), the governing unsteady 2D flow model is changed into two highly non-linear ODEs. After that, QP is applied to linearize the non-linear ODEs, and finally S2HWs and LIDQT are used to simulate the non-linear system of ODEs. In results and discussion section, the parameters Reynolds number (R), expansion ratio and Nusselt (Nu), and nanoparticle volume fractionφ are analysed with respect to velocity and temperature profiles. The proposed techniques are easy to implement for fluid and heat transfer (FHT) problems.


---

## Body

## 1. Introduction

In current century, the fluid flow and heat transfer (FFHT) are regarded as the extremely crucial problems of engineering and industries. This is the main reason that NFTT is a hot area of research amongst engineers and research community. It has wide applications in various scientific and industrial areas, such as power plant operations, manufacturing, lubrications, heat exchangers, micro-electromechanical technology, advanced nuclear energy systems, macromolecular science, and so on.Due to presence of nonlinear thermal radiation and viscous flow in the NFTT models, the analytical studies (in form of solutions) of these models are not an easy task. Secondly, the governing equations (GEs) of NFTT models are in the form of non-linear partial differential equations (PDEs) involving emerging parameters, which may be quite complex, and exact solutions of such type of problems are difficult to calculate. That is why a number of experimental and numerical studies on NFTT models have been conducted in the last few years [1–11].Carbon nanotubes (CNTs) with a cylindrical nanostructure were first studied by Venkataraman et al. and his associates [12] in 1991. Enough efforts were done to study the NFTT models through exact as well as classical numerical methods till date. Sheikholeslami Kandelousi [1] used the Runge–Kutta 4th order (RK4) scheme to study the characteristics of NFTT between two horizontal parallel plates. The study reflects that heat transmission (HT) growth increases as we increase Reynolds number when m=0, but the adverse tendency is inspected for different values of power law index m. Ahmed and his group [2] studied the magneto-hydrodynamic (MHD) fluid flow problem in a domain of rectangular shape. Basically, this technique is an iterative scheme and gives series solution, but this technique is quite complex due to series solution. Ahmed et al. [3] examined the NFTT in a rectangular channel with the help of the Galerkin method, while Hatami and his associates [4] studied the NFTT between two parallel plates by using the Galerkin and the least square methods. Khan and his group [5, 6] simulated the NFTT model with CNT-based nanofluids in non-parallel stretchable walls under the effects of velocity slip in a channel using differential transform (DT) and Runge–Kutta–Fehlberg schemes. Srinivas and his associates [13] analysed the thermal-diffusion and diffusion-thermo effects of a viscous fluid in a 2D channel between slowly expanding or contracting walls with weak permeability. Thus, there are many numerical techniques such as the similarity transformation (ST), homotopy analysis method (HAM), shooting method, Galerkin’s finite element method (FEM), Runge–Kutta method, Runge–Kutta–Fehlberg method, and so on [1–11] which have been proposed for these types of NFTT models. These types of models have been investigated in [14, 15]. For such kind of models, hybrid and new techniques are welcome; those are nicely implemented with less computation cost.The main contribution of this work is to cultivate two numerical algorithms based on S2HWs, LIDQT, and QP to simulate the unsteady nanofluid FHT in an asymmetric rectangular porous channel. S2HWs and LIDQT-based algorithms are powerful tools that fulfil all the required conditions of an efficient technique. Herein, the authors analysed the features of the unsteady NFTT model in a rectangular channel using S2HWs and LIDQT-based algorithms. In the simulation procedure, first of all, using similarity transformation (ST), the governing unsteady 2D flow model is transformed into two ODEs in which one is non-linear and other is linear. Then, QP is applied to linearize the non-linear ODEs, and finally S2HWs and LIDQT are used to simulate the non-linear system of ODEs. In results and discussion section, the parameters Reynolds number (R), expansion ratio and Nusselt (Nu), and nanoparticle volume fractionφ are analysed with respect to velocity and temperature profiles. The proposed techniques are easy to implement than the techniques available in the literature [3, 7, 16, 17] for fluid and heat transfer problems.


## 2. Description of the Problem and Mathematical Model

Consider a rectangular channel shown as in Figure1. The main end of the channel is supported by an impenetrable flexible membrane. Unsteady and laminar nanofluid is flowing through the channel. Other details can be found in [3]. For the geometrical explanation, take the middle part of the head end as original displayed in Figure 1. Here, we assume that the temperature T of the lower wall Al is greater than the upper wall Au. In these considerations, the basic mass, momentum, and energy governing equations of the problem are(1)∂U∂ξ+∂V∂η=0,(2)ρnf∂U∂t+U∂U∂ξ+V∂U∂η=−∂P∂ξ+μnf∂2U∂ξ2+∂2U∂η2,(3)ρnf∂V∂t+U∂V∂ξ+V∂V∂η=−∂P∂η+μnf∂2V∂ξ2+∂2V∂η2,(4)∂T∂t+U∂T∂ξ+V∂T∂η=knfρCpnf∂2T∂ξ2+∂2T∂η2.Figure 1
Geometry of the problem.Here,U is the horizontal component of velocity and V is the vertical component of velocity. The viscosity μnf, density ρnf, and effective thermal conductivity knf of nanofluids and effective heat capacity ρCpnf are in the following forms:(5)ρnf=1−ϕρf+ϕρCNT,μnf=μf1−ϕ2.5,ρCpnf=1−ϕρCpf+ϕρCpCNT,knfkf=1−ϕ+2ϕkCNT/kCNT−kflnkCNT+kf/2kf1−ϕ+2ϕkf/kCNT−kflnkCNT+kf/2kf.Equations (1)–(4) are connected with the BCs(6)iη=at:U=0,V=−Vu=−Aua˙,T=Tu,(7)iiη=−at:U=0,V=−Vl=−Ala˙,T=Tl,where Al and Au represent the penetrability of lower and top walls, respectively, and a˙=da/dt is the velocity of moving wall. The values of parameters of base fluid and nanoparticles are taken from Table 1.Table 1
Parameters of base fluid and nanoparticles.
ρkg/m3kW/mKCpj/kgKPrPure water997.10.61341796.2SWCNTs26006600425—MWCNTs16003000796—Now, to reduce the independent variables, choose the stream functionΦ=νξ/aGy,t, with y=η/at and kinematic viscosity ν=μnf/ρnf which satisfies continuity equation (1). Using these transformations, the velocity components become(8)U=∂Φ∂η=νξGyy,ta2,V=−∂Φ∂ξ=−VGy,ta,where Gy=∂G/∂y.Substitute the value of the velocity componentsU and V from (8) into (2) and (3) and eliminate the pressure term:(9)Gyyyy+A1αA2yGyyy+3Gyy+GGyyy−GyGyy−A1a2ν−1A2Gyyt=0.BCs given in equations (6) and (7) are reduced to(10)G−1,t=A1A2Rl,G1,t=A1A2Ru,Gy−1,t=0,Gy1,t=0,where A1/A2=ρnf/ρf/μnf/μf=1−ϕ+ϕρCNT/ρf1−ϕ2.5 are constant parameters, α=aa˙/ν is the wall expansion ratio which is very +ve and −ve for expansion and contraction, and Rl=aVl/ν and Ru=aVu/ν are the permeation RNs relevant to the bottom and top walls which are +ve due to injection of cooled fluid. Following the work of Uchida and Aoki and Pandit and Sharma [18, 19], to normalize equations (9) and (10), substitute g=G/R which leads to(11)g4+A1A2αyg3+3g2+Rgg3−g1g2=0,where g=gy and gr,r=1,2,3,4 denotes the derivatives of g with respect to y and the BCs are(12)g−1=A∗,g1=1,g1−1=0,g11=0.Now, to reduce equation (4) in one variable, use the following similarity transformation (ST):(13)T=Tu+Tl−Tuhy.Substitute equation (13) into energy equation (4) which leads the following ODE:(14)h2+A3A4PrRg+αyh1=0,with BCs(15)h−1=1,h1=0,where A3=1−ϕ+ϕρCpCNT/ρCpf, A4=knf/kf are constant parameters, and Pr is known as the Prandtl number. Finally, the non-dimensional Nusselt number (Nu) can be calculated as(16)Nul=−A4h−1,Nuu=−A4h1.


## 3. Quasilinearization Process for Linearization

To simulate equations (1)–(4), we reduce the equations into equations (11), (12), (14), and (15). Now, the next purpose is to solve non-linear ODE (11) and linear ODE (14). In this direction, the first step is to linearize non-linear ODE (11) by using QP proposed by Bellman and Kalaba [20]. The usefulness of the QP is that it linearizes the problem which is easy to solve than the non-linear problem. The non-linear terms of equation (11) are linearized as follows:(17)gg3n+1=gg3n+gn+1−gngn3+gn+13−gn3gn,orgg3n+1=gn+1gn3+gn+13gn−gg3n,(18)g1g2n+1=g1g2n+gn+11−gn1gn2+gn+12−gn2gn1,org1g2n+1=gn+11gn2+gn+12gn1−g1g2n,where n denotes the number of iterations of the quasilinearization technique.Remark 1.
The quasilinearization technique converges quadratically, i.e.,gn+1−gn≤Cgn−gn−12 where gn denotes the solution at nth step.
Now, use (17)and(18) in (11), and we have following linear system of ODE:(19)gn+14+A1A2αygn+13+3gn+12+A1A2Rgn+1gn3+gngn+13−gngn3−gn+11gn2−gn1gn+12+gn1gn2=0.


## 4. Wavelet Family and Uniform S2HWs

Wavelet-based techniques are powerful tools to solve differential equations (DEs). Generally, scaling functionχt and mother wavelet ψt are important to generate a family of wavelets. These two functions satisfy the following conditions:(20)∫−∞∞χ^ωωdω<∞,∫−∞∞ψ^ωωdω<∞,where χ^ω and ψ^ω are Fourier transforms of scaling function χt and mother wavelet ψt, respectively. The wavelet family can be generated by the following relation:(21)χkt=χt−k,ψj,kt=2j/2ψ2jt−k,where j,k are integers and these are called dilation and translation parameters, respectively.Any functiongt∈L2−∞,∞ can be represented into wavelets as(22)gt=∑kckχkt+∑k∑jdj,kψj,kt,where ck, dj,k are scaling and wavelet coefficients, respectively. These coefficients can be calculated by the orthogonal property of wavelets.


### 4.1. Uniform S2HWs

Uniform S2HWs are newly developed resources for the simulation of DEs and integral equations. These wavelets were proposed early by Alfred Haar in 1910. These wavelets are symmetric Daubechies wavelets of order 1 with vanishing 0 moments. Uniform S2HWs forx∈0,1 are defined as(23)ψjy=1,y∈ξ1j,ξ2j,−1,y∈ξ2j,ξ3j,0,elsewhere,where(24)ξ1j=km,ξ2j=k+0.5m,ξ3j=k+1m.In equation (24), m=2p,p=0,1,…,J; k=0,1,…,m−1 indicate the wavelet level and translation parameter, respectively. As wavelet level J increases, the wavelets become finer and finer. The relation between j,k,m is given as j=m+k+1. The scaling function for j=1 is defined in the following form:(25)ψ1y=χy=1,y∈0,1,0,otherwise.To solve ODEs, we need the integrals of uniform S2HWs (23)–(25). The βth order integrals are defined as follows:(26)pβ,jy=∫0y∫0tβ−1…∫0t1ψjtdtdt1…dtβ−1︸β−times=1α−1!∫0yy−tβ−1hitdt,β=1,2,…,n;j=1,2,…,2M.Solving integral (26) analytically, we have the following form:(27)pβ,jy=0,fory<ξ1j,y−ξ1jββ!,fory∈ξ1j,ξ2j,y−ξ1jβ−2y−ξ2jββ!,fory∈ξ2j,ξ3j,y−ξ1jβ−2y−ξ2jβ+y−ξ3jββ!,fory>ξ3j.


### 4.2. Uniform S2HWs for Simulation of the Problem

To simulate the system of ODEs (19) and (14), we approximate the highest order derivatives in the linear combination of uniform S2HWs as follows [21–27].(28)gn+14y=∑i=12Mdiψiy=d2MTψ2My,where d2MT=d1,d2,…,d2M represent unknown S2HW coefficients.Integrate successively equation (28) from 0 to y to obtain the lower-order derivatives. The lower-order derivatives in terms of S2HW operational matrices are as follows:(29)gn+13y=d2MTP1ψ2My+gn+130,gn+12y=d2MTP2ψ2My+ygn+130+gn+120,gn+11y=d2MTP3ψ2My+y22gn+130+ygn+120+gn+110,gn+1y=d2MTP4ψ2My+y36gn+130+y22gn+120+ygn+110+gn+10.Similarly, we can obtain the approximation in terms of S2HW operational matrices for equation (14) as follows:(30)h2y=∑i=12Meiψiy=e2MTψ2My,h1y=e2MTP1ψ2My+g10,hy=e2MTP2ψ2My+yg10+g0,where e2MT=e1,e2,…,e2M are unknown S2HW coefficients.First, apply the boundary conditions (12) and (15) in equations (29) and (30). After that, substitute the values of gn+14y,gn+13y,gn+12y,gn+11y,gn+1y from (29) in (19) and h2y,h1y,hy from (30) in (14) with yl=l−0.5/2M,l=1,2,…,2M, and we have the linear systems:(31a)Q1D=B1,(31b)Q2E=B2,where Q1, Q2 are known S2HW operational matrices of order 2M×2M and D, E are unknown S2HW coefficient matrices of order 2M×1.Solve the systems in (31a) and (31b) with the LU decomposition method for the S2HW coefficients dl,el,l=1,2,…,2M. Further, with the help of these S2HW coefficients and yl=l−0.5/2M,l=1,2,…,2M, we can obtain the approximate solutions (ASs) from (29) and (30) in the considered region.


### 4.3. Convergence Analysis of S2HW Approximation

This section describes the convergence criteria of the S2HWs via the following result.Lemma 1. (see [26]).
Letgy∈L2a,b and gy=∑l=1∞dlψly. If g′y≤L,l>0, then dl≤L2−3l−2/2.Theorem 1.
Letgy be the exact solution of (11) and G2My be simulated solution via S2HWs. Let ER=gy−G2My be the error; then,(32)ER≤C1M2−32J−11−2−3/2,where J is the level of resolution of S2HWs.Proof.
Givengy=∑l=1∞dlψly and G2My=∑l=12Mdlψly, then the error estimation at Jth level resolution is(33)ER2≤gy−G2My2=∑l=2M+1∞dlψly2=∑l=2M+1∞dlψly,∑l=2M+1∞dlψly.
Using the inner product, orthonormal condition of S2HWs, mean value theorem, and Lemma1, we have(34)ER2≤K2−32J−11−2−3/2∑l=2M+1∞dlKl.
AssumeC1=SuplKl, and again using Lemma 1, we conclude that(35)ER≤C1M2−32J−11−2−3/2.
This reflects that asJ increases, the error goes to zero. Hence, the solution via S2HWs converges to exact solution.


## 4.1. Uniform S2HWs

Uniform S2HWs are newly developed resources for the simulation of DEs and integral equations. These wavelets were proposed early by Alfred Haar in 1910. These wavelets are symmetric Daubechies wavelets of order 1 with vanishing 0 moments. Uniform S2HWs forx∈0,1 are defined as(23)ψjy=1,y∈ξ1j,ξ2j,−1,y∈ξ2j,ξ3j,0,elsewhere,where(24)ξ1j=km,ξ2j=k+0.5m,ξ3j=k+1m.In equation (24), m=2p,p=0,1,…,J; k=0,1,…,m−1 indicate the wavelet level and translation parameter, respectively. As wavelet level J increases, the wavelets become finer and finer. The relation between j,k,m is given as j=m+k+1. The scaling function for j=1 is defined in the following form:(25)ψ1y=χy=1,y∈0,1,0,otherwise.To solve ODEs, we need the integrals of uniform S2HWs (23)–(25). The βth order integrals are defined as follows:(26)pβ,jy=∫0y∫0tβ−1…∫0t1ψjtdtdt1…dtβ−1︸β−times=1α−1!∫0yy−tβ−1hitdt,β=1,2,…,n;j=1,2,…,2M.Solving integral (26) analytically, we have the following form:(27)pβ,jy=0,fory<ξ1j,y−ξ1jββ!,fory∈ξ1j,ξ2j,y−ξ1jβ−2y−ξ2jββ!,fory∈ξ2j,ξ3j,y−ξ1jβ−2y−ξ2jβ+y−ξ3jββ!,fory>ξ3j.


## 4.2. Uniform S2HWs for Simulation of the Problem

To simulate the system of ODEs (19) and (14), we approximate the highest order derivatives in the linear combination of uniform S2HWs as follows [21–27].(28)gn+14y=∑i=12Mdiψiy=d2MTψ2My,where d2MT=d1,d2,…,d2M represent unknown S2HW coefficients.Integrate successively equation (28) from 0 to y to obtain the lower-order derivatives. The lower-order derivatives in terms of S2HW operational matrices are as follows:(29)gn+13y=d2MTP1ψ2My+gn+130,gn+12y=d2MTP2ψ2My+ygn+130+gn+120,gn+11y=d2MTP3ψ2My+y22gn+130+ygn+120+gn+110,gn+1y=d2MTP4ψ2My+y36gn+130+y22gn+120+ygn+110+gn+10.Similarly, we can obtain the approximation in terms of S2HW operational matrices for equation (14) as follows:(30)h2y=∑i=12Meiψiy=e2MTψ2My,h1y=e2MTP1ψ2My+g10,hy=e2MTP2ψ2My+yg10+g0,where e2MT=e1,e2,…,e2M are unknown S2HW coefficients.First, apply the boundary conditions (12) and (15) in equations (29) and (30). After that, substitute the values of gn+14y,gn+13y,gn+12y,gn+11y,gn+1y from (29) in (19) and h2y,h1y,hy from (30) in (14) with yl=l−0.5/2M,l=1,2,…,2M, and we have the linear systems:(31a)Q1D=B1,(31b)Q2E=B2,where Q1, Q2 are known S2HW operational matrices of order 2M×2M and D, E are unknown S2HW coefficient matrices of order 2M×1.Solve the systems in (31a) and (31b) with the LU decomposition method for the S2HW coefficients dl,el,l=1,2,…,2M. Further, with the help of these S2HW coefficients and yl=l−0.5/2M,l=1,2,…,2M, we can obtain the approximate solutions (ASs) from (29) and (30) in the considered region.


## 4.3. Convergence Analysis of S2HW Approximation

This section describes the convergence criteria of the S2HWs via the following result.Lemma 1. (see [26]).
Letgy∈L2a,b and gy=∑l=1∞dlψly. If g′y≤L,l>0, then dl≤L2−3l−2/2.Theorem 1.
Letgy be the exact solution of (11) and G2My be simulated solution via S2HWs. Let ER=gy−G2My be the error; then,(32)ER≤C1M2−32J−11−2−3/2,where J is the level of resolution of S2HWs.Proof.
Givengy=∑l=1∞dlψly and G2My=∑l=12Mdlψly, then the error estimation at Jth level resolution is(33)ER2≤gy−G2My2=∑l=2M+1∞dlψly2=∑l=2M+1∞dlψly,∑l=2M+1∞dlψly.
Using the inner product, orthonormal condition of S2HWs, mean value theorem, and Lemma1, we have(34)ER2≤K2−32J−11−2−3/2∑l=2M+1∞dlKl.
AssumeC1=SuplKl, and again using Lemma 1, we conclude that(35)ER≤C1M2−32J−11−2−3/2.
This reflects that asJ increases, the error goes to zero. Hence, the solution via S2HWs converges to exact solution.


## 5. LIDQT for Simulation of the Problem

In this section, first we will give brief details of LIDQT and then use the technique for the simulation purpose.


### 5.1. Overview of LIDQT

The differential quadrature method (DQM) was introduced by Bellman et al. [28] in 1972, and it is based on Lagrange’s interpolation. That is why it is called LIDQT. LIDQT became famous for solving differential equations (DEs) after 1972. The main use of LIDQT is to approximate the derivatives of unknown functions in the problems identical to integral quadrature rules. To understand the idea of the LIDQT, consider a sufficiently smooth function fy defined on the domain α,β and let the interval be divided into the grid nodes as α=y1,<y2,…,<yp=β. By using the DQ technique, the qth order derivative fqy is approximated at y=yi in the following way:(36)fqyi≅∑j=1pωijqfyj,i=1,2,…,p,where p is the number of nodes in the interval α,β, yi is an arbitrary node, and ωijq are the weighting coefficients (WCs) of qth order derivative. In DQ technique methodology, the first step is to calculate the WCs ωijq. In the literature, there are many ways [16, 17, 20, 28–33] to compute the WCs. Here, we used the following Lagrange’s test functions (LTFs) to compute the WCs:(37)sly=Nyy−ykN1yk,l=1,2,…,p,where Ny=y−y1y−y2…y−yp, N1yk=∏j=1j≠kpyk−yj.Substitute (37) into (36), and we have(38)slqyi≅∑j=1pωijqsyj,i=1,2,…,p.Using the methodology used by Chang [34], we have the following 1st and higher-order WCs:(39)ωlk1=N1ylyl−ykN1yk,forl≠k,ωkk1=−∑l=1,l≠kpalk1,k=l,l,k=1,2,…,p,(40)ωlkq=qωlk1ωllq−1−ωlkq−1yl−yk,forl≠k,ωkkq=−∑l=1,l≠kpωlkq,k=l,l,k=1,2,…,p.


### 5.2. Simulation of the Problem with LIDQT

Now, we will apply the above-explained LIDQT on equations (19) and (14) for numerical simulation with the BCs (12) and (15), respectively. Here, we choose the non-uniform Chebyshev–Gauss–Lobatto (CGL) nodes yk=1/21−cosk−1π/p−1L,k=1,2,…,p for discretization, where L is length of the domain. Discretize the derivatives at the nodes yi,i=1,2,…,p in equations (19) and (14), and we have(41)∑k=1pωik4fk+A1A2αyi∑k=1pωik3fk+3∑k=1pωik2fk+Rfif0i3+f0i∑k=1pωik3fk−f0i2∑k=1pωik1fk−f0i1∑k=1pωik1fk=A1A2Rf0i∑k=1pωik3f0k−f0i2∑k=1pωik1f0k,(42)∑k=1pωik2hk+A3A4prRfi+αyi∑k=1pωik1hk=0.Here, for simplicity, we assumegn+1=f,gn=f0 in equation (41). Equations (41) and (42) are two linear systems of algebraic equations, which are solved after implementation of BCs by the LU decomposition method simultaneously.


### 5.3. Stability Criteria of LIDQT

Stability analysis is a major issue for solving the fluid and heat problems. The stability of the LIDQT depends on the eigenvalues (EVs) of WC matrices, and these EVs totally depend on the distribution of node points used for domain discretization. This concept had been studied and shown by Chang [34]. In this book, Shu et al. showed that the uniform node distribution does not produce stable solution, but the stable solution can be acquired by using CGL nodes which we have used for simulation purpose.


## 5.1. Overview of LIDQT

The differential quadrature method (DQM) was introduced by Bellman et al. [28] in 1972, and it is based on Lagrange’s interpolation. That is why it is called LIDQT. LIDQT became famous for solving differential equations (DEs) after 1972. The main use of LIDQT is to approximate the derivatives of unknown functions in the problems identical to integral quadrature rules. To understand the idea of the LIDQT, consider a sufficiently smooth function fy defined on the domain α,β and let the interval be divided into the grid nodes as α=y1,<y2,…,<yp=β. By using the DQ technique, the qth order derivative fqy is approximated at y=yi in the following way:(36)fqyi≅∑j=1pωijqfyj,i=1,2,…,p,where p is the number of nodes in the interval α,β, yi is an arbitrary node, and ωijq are the weighting coefficients (WCs) of qth order derivative. In DQ technique methodology, the first step is to calculate the WCs ωijq. In the literature, there are many ways [16, 17, 20, 28–33] to compute the WCs. Here, we used the following Lagrange’s test functions (LTFs) to compute the WCs:(37)sly=Nyy−ykN1yk,l=1,2,…,p,where Ny=y−y1y−y2…y−yp, N1yk=∏j=1j≠kpyk−yj.Substitute (37) into (36), and we have(38)slqyi≅∑j=1pωijqsyj,i=1,2,…,p.Using the methodology used by Chang [34], we have the following 1st and higher-order WCs:(39)ωlk1=N1ylyl−ykN1yk,forl≠k,ωkk1=−∑l=1,l≠kpalk1,k=l,l,k=1,2,…,p,(40)ωlkq=qωlk1ωllq−1−ωlkq−1yl−yk,forl≠k,ωkkq=−∑l=1,l≠kpωlkq,k=l,l,k=1,2,…,p.


## 5.2. Simulation of the Problem with LIDQT

Now, we will apply the above-explained LIDQT on equations (19) and (14) for numerical simulation with the BCs (12) and (15), respectively. Here, we choose the non-uniform Chebyshev–Gauss–Lobatto (CGL) nodes yk=1/21−cosk−1π/p−1L,k=1,2,…,p for discretization, where L is length of the domain. Discretize the derivatives at the nodes yi,i=1,2,…,p in equations (19) and (14), and we have(41)∑k=1pωik4fk+A1A2αyi∑k=1pωik3fk+3∑k=1pωik2fk+Rfif0i3+f0i∑k=1pωik3fk−f0i2∑k=1pωik1fk−f0i1∑k=1pωik1fk=A1A2Rf0i∑k=1pωik3f0k−f0i2∑k=1pωik1f0k,(42)∑k=1pωik2hk+A3A4prRfi+αyi∑k=1pωik1hk=0.Here, for simplicity, we assumegn+1=f,gn=f0 in equation (41). Equations (41) and (42) are two linear systems of algebraic equations, which are solved after implementation of BCs by the LU decomposition method simultaneously.


## 5.3. Stability Criteria of LIDQT

Stability analysis is a major issue for solving the fluid and heat problems. The stability of the LIDQT depends on the eigenvalues (EVs) of WC matrices, and these EVs totally depend on the distribution of node points used for domain discretization. This concept had been studied and shown by Chang [34]. In this book, Shu et al. showed that the uniform node distribution does not produce stable solution, but the stable solution can be acquired by using CGL nodes which we have used for simulation purpose.


## 6. Simulation and Remarks

The main goal of this segment is to present the simulated results of model via the developed numerical algorithms based on S2HWs, LIDQT, and QP. The computational work is done with 31 and 64 grid points for LIDQT and S2HWs, respectively. The maximum absolute errorL∞ and relative error (RE) are computed by the formulas(43)L∞=uanalytic−uappx.∞=maxjujanalytic−ujappx.,RE=Δx∑k=1Nukanalytic−ukappx.2,where uanalytic and uappx. denote analytical and approximated solutions, respectively. The simulation part has been done by writing MATLAB routines.First, to accredit the developed algorithms for the model (11), let us assume gy=y4 as exact solution (ES) of equation (11). This assumption coverts homogeneous equation (11) into a non-homogeneous equation as follows:(44)g4+A1A2αyg3+3g2+Rgg2−g1g2=χy,where(45)χy=24+60αA1A2y2−24RA1A2y5,with BCs(46)g−1=1,g1=1,g1−1=−4,g11=4.The main reason of this assumption is to validate the developed algorithms for the problem since now we can compute absolute (L∞) and relative errors (REs) of the problem. After that, systems (44)–(46) are solved by the developed algorithms. Obtained outcomes are given in Table 2 and Figures 2 and 3. Table 2 makes a comparison of S2HWs and LIDQT in form of the errors L∞ and REs with parameters α=−1,R=−1.5,φ=0.9, and the table concludes that both algorithms give good accuracy and efficiency. The computational costs of the algorithms S2HWs and LIDQT are 1.325913 s and 1.23521 s, respectively. Figures 2 and 3 compare the exact, S2HW, and LIDQT solutions, and these figures reflect that S2HW and LIDQT solutions are almost the same. Thus, the table and figures justified the competence and good convergence rate of the developed algorithms. Figures 2 show the exact, S2HW, and LIDQT solutions for fixed values of EC ratio α=−1, Reynolds number R=−1.5, and NVF φ=0.9. The behaviour of gradient of velocity (GV) with respect to to y for φ=0.9,α=1, and R=1 is shown in Figure 3. Figures 4 and 5 show the parameters R=0.1,φ=0.1 with the BCs A∗=−0.2,0.2, respectively. Figure 4 reflects the growth or decay of velocity with different EC ratio α while Figure 5 depicts velocity profile with respect to different RNR. Figures 2–4 are plotted via S2HWs and LIDQT, and we conclude that results are almost same with S2HWs and LIDQT algorithms. That is why after Figure 3, we have reported results via S2HW algorithm only. Figures 6 and 7 portray the variation of GV behaviour g1y, for R=1 and R=−1 with EC ratio between α∈−2.5,2.5 for SWCNTs notified by concrete lines and MWCNTs signified by dashed lines. Figures 6 and 7 conclude that the maximum velocity lies at the centre of the channel, while near the wall, GV profile decreases as EC ratio α increases. Figure 8 discusses the case when contraction is combined with injection and suction, and the figure demonstrates that the velocity decreases with increasing value of R. Figure 9 reveals that as we decrease the value of R, the profile of GV goes down till R=1 and again goes up for R=0 in both the cases of SWCNTs and MWCNTs. Figures 2–9 show the impact of parameters on velocity gy and GV g′y profiles. Figure 10 shows the effect of NVF φ on temperature profile hy of problem (14). It is concluded that enhancement in NVF φ causes the rapid increase of temperature for MWCNTs as compared to SWCNTs. Figure 11 demonstrates that both SWCNTs and MWCNTs exhibit the same phenomenon when altering the value of Reynolds number for lower and upper walls in problems (11)–(14).Table 2
Absolute (L∞) and relative errors (REs) of problem (45) via S2HWs and LIDQT for α=−1,R=−1.5,φ=0.9.
yS2HWsLIDQTL∞REL∞RE−0.13.85E − 183.85E − 142.75E − 175.02E − 14−0.31.73E − 182.14E − 167.63E − 174.28E − 16−0.51.21E − 162.23E − 162.98E − 168.12E − 16−0.70.0000.0006.38E − 160.000−0.94.15E  −175.62E − 161.11E − 150.0000.18.58E − 188.57E − 145.74E − 177.74E − 140.31.56E − 171.93E − 153.99E − 171.50E − 150.52.78E − 164.44E − 160.0003.33E − 160.72.78E − 161.16E − 165.55E − 171.16E − 160.91.11E − 161.69E − 160.0000.000Figure 2
Exact, S2HW, and LIDQT solutions of problem (44) for the parameters φ=0.9,R=−1.5,α=−1 (a comparison).Figure 3
Exact, S2HW, and LIDQT solutions of problem (44) for the parameters R=1.5,α=1,φ=0.9 (a comparison).Figure 4
Effect of expansion ratioα over the profile g (y) via (a) S2HWs and (b) LIDQT for the parameters A∗=0.2,R=0.1,φ=0.1.
(a)(b)Figure 5
Effect of Reynolds numberR over the profile g (y) via S2HWs for the parameters (a) A∗=−0.2,α=−2,φ=0.1 and (b) A∗=0.2,α=−2,φ=0.1.
(a)(b)Figure 6
Impact of EC ratioα over the profile g1y for SWCNTs (concrete lines) and MWCNTs (dashed lines) via S2HWs for the parameters A∗=0.2,φ=0.1,R=1.Figure 7
Impact of RC ratioα over profile g1y for SWCNTs (dashed lines) and MWCNTs (concrete lines) via S2HWs for the parameters A∗=0.2,φ=0.1,R=−1.Figure 8
Impact ofR over the profileg1y for SWCNTs (concrete lines) and MWCNTs (dashed lines) via S2HWs for the parameters A∗=0.2,φ=0.1,α=−2.5.Figure 9
Impact ofR over the profileg1y for SWCNTs (concrete lines) and MWCNTs (dashed lines) via S2HWs for the parameters A∗=0.2,φ=0.1,α=3.Figure 10
Impact on temperature profilehy via S2HWs for the parameters A∗=−0.2,R=2.5,α=−1.5,Pr=6.2.Figure 11
Nu vs.φ impact of R via S2HWs for the parameters α=−2.5,A∗=−0.2,Pr=6.2.


## 7. Concluding Remarks

This work analysed the nanofluid FHT model by developing two numerical algorithms based on QP, S2HWs, and LIDQT. The developed algorithms are influenced by study done in [21, 22, 29] for fluid problems. Further, the impact of NVF, Reynolds number, and EC ratio has been investigated on the nanofluid FHT. The precision of developed algorithms has been tested by choosing the ES of equation (11) and by finding the L∞ and REs. It is concluded that SWCNTs have an exalted rate of heat transmission as compared to MWCNTs. Enhancement in the NVF does not influence the velocity outline, but it has plausible impact on temperature profile. Further, both single as well as multi-walled carbon nanotubes present the same demeanour for different values of Reynolds number for lower and upper walls. The proposed algorithms are easy to implement than the techniques used in [3, 7, 16, 17] for fluid and heat transfer problems. The proposed algorithms with some modifications can be useful for different types of fluid flow problems and higher-order non-linear mathematical models.



---

*Source: 5015018-2022-03-15.xml*