# Coiflet Wavelet-Homotopy Solution of Channel Flow due to Orthogonally Moving Porous Walls Governed by the Navier–Stokes Equations

Authors: Qing-Bo Chen; Hang Xu
Journal: Journal of Mathematics
 (2020)

---

## Abstract

A newly computational method based on the Coiflet wavelet and homotopy analysis method is developed, which inherits the great nonlinear treatment of the homotopy analysis technique and the local high-precision capability of the wavelet approach, to give solutions to the classic problem of channel flow with moving walls. The basic principle of this suggested technique and the specific solving process are presented in detail. Its validity and efficiency are then checked via rigid comparisons with other computational approaches. It is found that the homotopy-based convergence-control parameter and the wavelet-based resolution level of Coiflet are two effective ways to improve on accuracies of solutions.


---

## Body

## 1. Introduction

Liquid flow in porous channels has received extensive attention owing to its wide industrial and engineering applications. The earliest study on channel flow with porous walls was performed by Berman [1], who obtained a perturbation solution for small Reynolds number. Morduchow [2] formulated a closed-form solution of Berman’s problem [1] over the full injection range by means of the method of averages. Berman’s solution was then extended by Proudman [3] to the large Reynolds number case and further generalized by Wah [4] for both small and large Reynolds number cases. Other extended studies have been done by various researchers on different aspects such as the theoretical solution at high suction Reynolds number [5], experimental observation for injection and suction effects [6], exact solution at high suction or injection velocity [7], complete solution for large injection with the viscous layer effect [8], and high-precision series solution for the uniformly porous tube [9]. A significant extension was given by Brady and Acrivos [10] for inclusion of the accelerating walls into the classic channel flow problem. Watson et al. [11] further investigated the case that the channel walls undergo both acceleration and porosity. Dauenhauer and Majdalani [12] examined the unsteady channel flow caused by uniformly time-dependent wall expansion ratio in the presence of porous and moving walls. Their work was further extended by Xu et al. [13] to the nonuniform time-dependent wall expansion ratio case.In fact, most of fluid flow problems are nonlinear problems, which bring great difficulties to give their solutions. Various analytical or numerical methods have been proposed to solve these nonlinear flow problems such as the Runge–Kutta method [14, 15], the Akbari–Ganji method (AGM) [16, 17], the homotopy perturbation method (HPM), the collocation method (CM), the finite element method (FEM) [18], the Galerkin finite element method (GFEM) [19], and other numerical methods [20, 21]. Homotopy method is widely used because of its superiority in solving nonlinear problems. At the same time, many homotopy-based methods have also been proposed. Liu studied a parameter identification problem of two-phase porous media fractional flow equations [22] and then solved a nonlinear convection-diffusion equation [23] by a wavelet multiscale-homotopy method. In addition, based on the homotopy method, a multigrid-homotopy method was used to solve a nonlinear inverse problem governed by a convection-diffusion equation [24], and the numerical identification of diffusion parameters in a nonlinear convection-diffusion equation was investigated by an adaptive homotopy perturbation method [25].Recently, the homotopy analysis method (HAM) [26] was developed for handling strong nonlinear problems, which possesses several advantages such as freedom to choose solution expressions, independence of small or large parameters, and adjustment of solution convergence. It has been used to solve many kinds of nonlinear problems in different fields [27, 28]. However, solutions by this technique are purely based on analytical computation; as a result, it is very time consuming when nonlinear problems are very complicated. On the contrary, the Coiflet wavelet technique [29] received applausive attention owing to its excellent features including feasible wavelet transformation, high computational efficiency, and superior expression of frequency and time information. This technique was initially used in signal analysis [30] and recently extended by Zhou et al. [31, 32] and Liu et al. [33, 34] to some classic nonlinear homogeneous problems in applied mathematics and mechanics. Very recently, the homotopy-based Coiflet wavelet method was developed by Yang and Liao [35, 36] for capturing solutions with high precision to nonlinear equations. One-dimensional and two-dimensional Bratus equations subjected to homogeneous boundary conditions were considered in their work as illustrative examples. Their proposed technique was then extended by Yu and his collaborators to solve the nonlinear plate bending problem [37] and cavity flow problem [38] with nonhomogeneous boundary conditions.Though several successful applications of the homotopy-Coiflet wavelet technique have been achieved, many aspects of this approach are not yet clear. Further efforts are still needed to improve on its computational capability and efficiency. The aim of this paper is to extend the homotopy-Coiflet wavelet method to solve the classic porous channel flow problem governed by a fourth-order strong nonlinear system subjected to nonhomogeneous boundary conditions. The basic idea of this approach will be exhibited in detail. A range of wall expansion ratios for both injection and suction cases, especially for large injection and suction, will be taken into account. The accuracy and efficiency of our results will be then checked via a rigid comparison with previous studies.


## 2. Mathematical Description

As shown in Figure1, liquid flow is driven through a two-dimensional rectangular channel bounded by a porous wall at y=±a. The channel walls undergo uniform expansion or contraction. The constant ratio of fluid injection and suction through the two opposing channel walls is applied. The fixed Cartesian coordinates x,y in the space are used, with the x-axis extending along the length of the channel and the y-axis in the wall-normal direction.Figure 1
Physical sketch.The fourth-order nonlinear system describing fluid flow in this porous channel with moving walls is written as (refer to Dauenhauer and Majdalani [12])(1)F″″+FF‴+3α−F′F″+αηF‴=0,subject to the boundary conditions(2)F0=0,F″0=0,F1=R,F′1=0,where prime denotes the derivation with η η=y/a, F is the dimensionless stream function and the function of η, α is the wall expansion and extraction ratio, and R is the Reynolds number with positive R representing the injection case and vice versa. Note that this equation is derived directly from the Navier–Stokes equation via scaling and vorticity transformations without reduction. As a result, its solution can be deemed as the exact one of this famous equation.


### 2.1. Linearization

Here, the homotopy analysis method is employed to linearize nonlinear equation (1) and its boundary conditions (2). In doing so, we construct the following zeroth-order deformation equation:(3)1−qℒϕη;q−F0η=c0qNϕη;q,subject to the boundary conditions(4)ϕ0;q=0,∂2ϕη;q∂η2η=0=0,ϕ1;q=R,∂ϕη;q∂ηη=1=0,where q is an embedding parameter evolving from 0 to 1, ℒ is the auxiliary linear operator, F0η is the initial guess that satisfies all given boundary conditions, c0 is the convergence-control parameter, ϕη;q is the mapping function of Fη, and Nϕη;q is the nonlinear operator, which is defined by(5)Nϕη;q=∂4ϕη;q∂η4+ϕη;q∂3ϕη;q∂η3+3α∂2ϕη;q∂η2−∂ϕη;q∂η∂2ϕη;q∂η2+αη∂3ϕη;q∂η3.Since the computational domain is0,1, the highest-order derivative of equation (1) is fourth-order, and the quantity of boundary conditions is four; we choose ℒϕx;q=∂4ϕx;q/∂x4 as the auxiliary linear operator for this problem.From equation (3), it is readily known ϕη;0=F0η and ϕη;1=Fη. To get the exact solution Fη, we expand the solution ϕη;q of equation (3) by(6)ϕη;q=F0η+∑m=1+∞Fmηqm,where Fmη is the high-order subfunction. Clearly, Fη=ϕη;1 is acquired at q=1.Note that the initial guess should satisfy the given boundary conditions. According to equation (4), the initial guess F0η is chosen as(7)F0η=12Rη3−η2.To determineFmη, we construct the high-order deformation equation, by differentiating equation (3) m times with respect to q, then dividing by m!, and finally setting q=0, as(8)ℒFmη−χmFm−1η=c0RmF⟶m−1η,where(9)χm=0,m=1,1,m>1Fmη=1m!∂mϕη;q∂qmq=0,RmF⟶m−1η=1m−1!∂m−1Nϕη;q∂qm−1q=0,in which F⟶m−1η denotes a group of subfunctions F0η,F1η,F2η,…,Fm−1η.The corresponding high-order boundary conditions then take the form(10)Fm0=0,Fm″0=0,Fm1=0,Fm′1=0.


### 2.2. Coiflet Wavelet Expression

Different from the pure HAM technique, we select the Coiflet wavelet function rather than the power series given by Xu et al. [13] as the solution expression of Fη.Following Wang [39] and Liu [40], an arbitrary function fx is approximated by Coiflet wavelet as(11)fx≈∑k=2−3N+M12j−1+M1fk2jφ2jx−k+M1,where N is the vanishing moment, M1 is the first-order vanishing moment, and φ is the wavelet basis.In practical applications of wavelet approximation in boundary value problems, boundary modification is needed for handling boundary conditions. The modified boundary technique [32, 34, 39] for Coiflet wavelet is constructed as(12)φj,kx=∑i=2−3N+M1−1T0,ki2jφj,i+φj,k,k∈0,3,φj,k,k∈4,2j−4,∑i=2j+12j−1+M1T1,2j−ki2jφj,i+φj,k,k∈2j−3,2j,in which(13)φj,i=ϕ2jx−i+M1,φj,k=ϕ2jx−k+M1,T0,kx=∑i=03p0,i,ki!xi,T1,kx=∑i=03p1,i,ki!x−1i.Note thatp0,i,k and p1,i,k are coefficients determined by the following coefficient matrices P0 and P1(14)P0=1000−1163−32132−54−1−13−31,P1=1000116−332−132−54−11−33−1,via the relation P0=2−ijp0,i,k and P1=2−ijp1,i,k with i=0,1,2,3. Taking the boundary conditions (4) into consideration, it is known that Fη satisfies the inhomogeneous Dirichlet and homogeneous Neumann types of boundary conditions. Hence, following Yu et al. [37, 38, 41–43], by modifying the coefficient matrices, the Coiflet wavelet basis is selected as(15)hj,kη=φj,kηp0,2,k⟶0,p1,1,k⟶0.Furthermore,Fmη and Rmη can be expressed by Coiflet wavelet as(16)Fmη≈∑k=02jFmk2jhj,kη=∑k=12j−1Fmk2jhj,kηRmF⟶m−1η≈∑k=02jRmk2jφj,kη,where(17)Rmk2j=Fm−1″″k2j+3αFm−1″k2j+αηFm−1‴k2j+∑s=0m−1Fsk2jFm−1−s‴k2j−Fs′k2jFm−1−s″k2j.Note that the boundary conditions on both sides of the Coiflet wavelet expansion ofFmη are rectified due to the boundary condition constraints at η=0 and η=1.Generally, from Coiflet wavelet approximation (16), the nth-order derivatives with respect to η are obtained by the following formula as(18)dnFmηdηn≈∑k=12j−1Fmk2jdnhj,kηdηn.Therefore,mth-order deformation equation (8) is expressed as(19)∑k=12j−1Fmk2j−χmFm−1k2jℒhj,kη=c0∑k=02jRmk2jφj,kη.Following Yang and Liao [35, 36], the Galerkin method is employed to solve mth-order deformation equation (19). Multiplying both sides of equation (19) by hj,lηl=1,2,3,…,2j−1 and then integrating on the domain 0,1, we obtain the mth-order wavelet-Galerkin deformation equation as(20)ATF^m−χmF^m−1=c0BTR^m,where(21)F^m=Fm12j,Fm22j,…,Fm2j−12jT,R^m=Rm02j,Rm12j,…,Rm2j2jT,and the connection coefficient matrices are defined as(22)A=ak,l=Γk,1j,4k,l=1k,l=2j−1,B=bk,l=Γ¯k,1j,0k=0,l=1k=2j,l=2j−1,in which Γk,lj,n and Γ¯k,lj,n are the generalized connection coefficients defined by(23)Γk,1j,n=∫01dnhj,kξdxnhj,lξdξ,Γ¯k,1j,n=∫01dnφj,kξdxnhj,lξdξ.Note that the generalized connection coefficients change as the auxiliary linear operatorℒ and boundary conditions vary.By solving themth-order wavelet-Galerlin equations from the order m=1,2,…, the wavelet homotopy solution can be completely determined.


## 2.1. Linearization

Here, the homotopy analysis method is employed to linearize nonlinear equation (1) and its boundary conditions (2). In doing so, we construct the following zeroth-order deformation equation:(3)1−qℒϕη;q−F0η=c0qNϕη;q,subject to the boundary conditions(4)ϕ0;q=0,∂2ϕη;q∂η2η=0=0,ϕ1;q=R,∂ϕη;q∂ηη=1=0,where q is an embedding parameter evolving from 0 to 1, ℒ is the auxiliary linear operator, F0η is the initial guess that satisfies all given boundary conditions, c0 is the convergence-control parameter, ϕη;q is the mapping function of Fη, and Nϕη;q is the nonlinear operator, which is defined by(5)Nϕη;q=∂4ϕη;q∂η4+ϕη;q∂3ϕη;q∂η3+3α∂2ϕη;q∂η2−∂ϕη;q∂η∂2ϕη;q∂η2+αη∂3ϕη;q∂η3.Since the computational domain is0,1, the highest-order derivative of equation (1) is fourth-order, and the quantity of boundary conditions is four; we choose ℒϕx;q=∂4ϕx;q/∂x4 as the auxiliary linear operator for this problem.From equation (3), it is readily known ϕη;0=F0η and ϕη;1=Fη. To get the exact solution Fη, we expand the solution ϕη;q of equation (3) by(6)ϕη;q=F0η+∑m=1+∞Fmηqm,where Fmη is the high-order subfunction. Clearly, Fη=ϕη;1 is acquired at q=1.Note that the initial guess should satisfy the given boundary conditions. According to equation (4), the initial guess F0η is chosen as(7)F0η=12Rη3−η2.To determineFmη, we construct the high-order deformation equation, by differentiating equation (3) m times with respect to q, then dividing by m!, and finally setting q=0, as(8)ℒFmη−χmFm−1η=c0RmF⟶m−1η,where(9)χm=0,m=1,1,m>1Fmη=1m!∂mϕη;q∂qmq=0,RmF⟶m−1η=1m−1!∂m−1Nϕη;q∂qm−1q=0,in which F⟶m−1η denotes a group of subfunctions F0η,F1η,F2η,…,Fm−1η.The corresponding high-order boundary conditions then take the form(10)Fm0=0,Fm″0=0,Fm1=0,Fm′1=0.


## 2.2. Coiflet Wavelet Expression

Different from the pure HAM technique, we select the Coiflet wavelet function rather than the power series given by Xu et al. [13] as the solution expression of Fη.Following Wang [39] and Liu [40], an arbitrary function fx is approximated by Coiflet wavelet as(11)fx≈∑k=2−3N+M12j−1+M1fk2jφ2jx−k+M1,where N is the vanishing moment, M1 is the first-order vanishing moment, and φ is the wavelet basis.In practical applications of wavelet approximation in boundary value problems, boundary modification is needed for handling boundary conditions. The modified boundary technique [32, 34, 39] for Coiflet wavelet is constructed as(12)φj,kx=∑i=2−3N+M1−1T0,ki2jφj,i+φj,k,k∈0,3,φj,k,k∈4,2j−4,∑i=2j+12j−1+M1T1,2j−ki2jφj,i+φj,k,k∈2j−3,2j,in which(13)φj,i=ϕ2jx−i+M1,φj,k=ϕ2jx−k+M1,T0,kx=∑i=03p0,i,ki!xi,T1,kx=∑i=03p1,i,ki!x−1i.Note thatp0,i,k and p1,i,k are coefficients determined by the following coefficient matrices P0 and P1(14)P0=1000−1163−32132−54−1−13−31,P1=1000116−332−132−54−11−33−1,via the relation P0=2−ijp0,i,k and P1=2−ijp1,i,k with i=0,1,2,3. Taking the boundary conditions (4) into consideration, it is known that Fη satisfies the inhomogeneous Dirichlet and homogeneous Neumann types of boundary conditions. Hence, following Yu et al. [37, 38, 41–43], by modifying the coefficient matrices, the Coiflet wavelet basis is selected as(15)hj,kη=φj,kηp0,2,k⟶0,p1,1,k⟶0.Furthermore,Fmη and Rmη can be expressed by Coiflet wavelet as(16)Fmη≈∑k=02jFmk2jhj,kη=∑k=12j−1Fmk2jhj,kηRmF⟶m−1η≈∑k=02jRmk2jφj,kη,where(17)Rmk2j=Fm−1″″k2j+3αFm−1″k2j+αηFm−1‴k2j+∑s=0m−1Fsk2jFm−1−s‴k2j−Fs′k2jFm−1−s″k2j.Note that the boundary conditions on both sides of the Coiflet wavelet expansion ofFmη are rectified due to the boundary condition constraints at η=0 and η=1.Generally, from Coiflet wavelet approximation (16), the nth-order derivatives with respect to η are obtained by the following formula as(18)dnFmηdηn≈∑k=12j−1Fmk2jdnhj,kηdηn.Therefore,mth-order deformation equation (8) is expressed as(19)∑k=12j−1Fmk2j−χmFm−1k2jℒhj,kη=c0∑k=02jRmk2jφj,kη.Following Yang and Liao [35, 36], the Galerkin method is employed to solve mth-order deformation equation (19). Multiplying both sides of equation (19) by hj,lηl=1,2,3,…,2j−1 and then integrating on the domain 0,1, we obtain the mth-order wavelet-Galerkin deformation equation as(20)ATF^m−χmF^m−1=c0BTR^m,where(21)F^m=Fm12j,Fm22j,…,Fm2j−12jT,R^m=Rm02j,Rm12j,…,Rm2j2jT,and the connection coefficient matrices are defined as(22)A=ak,l=Γk,1j,4k,l=1k,l=2j−1,B=bk,l=Γ¯k,1j,0k=0,l=1k=2j,l=2j−1,in which Γk,lj,n and Γ¯k,lj,n are the generalized connection coefficients defined by(23)Γk,1j,n=∫01dnhj,kξdxnhj,lξdξ,Γ¯k,1j,n=∫01dnφj,kξdxnhj,lξdξ.Note that the generalized connection coefficients change as the auxiliary linear operatorℒ and boundary conditions vary.By solving themth-order wavelet-Galerlin equations from the order m=1,2,…, the wavelet homotopy solution can be completely determined.


## 3. Validation

### 3.1. Validation in the PDE Case

In this section, we will present a 2D linear problem to verify the accuracy and efficiency of our proposed method. The governing equation considered here is as follows:(24)∂2f∂x2+∂2f∂y2+2π2f1=0,where(25)f1x,y=cosπxcosπy.The corresponding homogeneous Neumann boundary conditions are(26)∂f∂xx=0=0,∂f∂xx=1=0,∂f∂yy=0=0,∂f∂yy=1=0.Note that equation (24) has the analytical solution(27)fex,y=cosπxcosπy.In our approach, by modifying the coefficient matrices, the Coiflet wavelet basis is expressed by(28)hj,kx=φj,kxp0,1,i⟶0,p1,1,i⟶0,hj,ly=φj,lyp0,1,i⟶0,p1,1,i⟶0.Furthermore,fx,y can be expressed by Coiflet wavelet as(29)fx,y≈∑k=02j∑l=02jfk2j,l2jhj,kxhj,ly,f1x,y≈∑k=02j∑l=02jfk2j,l2jφj,kxφj,ly.Therefore, equation (24) is expressed as(30)∑k=02j∑l=02jfk2j,l2jd2hj,kxdx2hj,ly+∑k=02j∑l=02jfk2j,l2jhj,kxd2hj,lydy2+2π2∑k=02j∑l=02jf1k2j,l2jφj,kxφj,ly=0.As described above, multiplying both sides of equation (30) by hj,nxhj,myn,m=0∼2j and then applying the multiple integrals on the domain 0,12, we obtain the wavelet-Galerkin equation as(31)A1˜T+A2˜T⋅vecf+B˜T⋅vecf1=0,where vec(⋅) is a straighten vector operator which pulls a matrix into a vector row by row and f,f1 is a 2j+1×2j+1 matrix defined by(32)f=fk,l=fk2j,l2j,f1=fk,l1=f1k2j,l2j,where the coefficient matrices are obtained through(33)A1˜=Γk,nj,2⊗Γl,mj,0,A2˜=Γk,nj,0⊗Γl,mj,2,B˜=2π2⋅Γ¯k,nj,0⊗Γ¯l,mj,0.Here,⊗ denotes the operator of the Kronecker tensor product, A1˜,A2˜,andB˜ are 2j+1×2j+1 matrices, and the generalized connection coefficients Γ and Γ¯ are defined as equation (23).To illustrate the error level of our solution, we define the average square errorErrSQ based on the exact solution (27) as follows:(34)ErrSQ=12j+12∑k=02j∑l=02jfk2j,l2j−fek2j,l2j2.The solution of equation (24) is shown in Figure 2, whose residual error is listed in Table 1. It is found that the average square error decreases rapidly with the resolution level j increasing. For j=7, especially, the error reaches to 10−9 within 1 minute CPU time. This completely demonstrates the accuracy and great efficiency of our proposed method.Figure 2
The wavelet solution of our example problem withj=6.Table 1
The computational time and average square error ErrSQ for equation (24).
Resolution levelj=4j=5j=6j=7CPU time (sec.)0.43490.63841.715332.3233ErrSQ2754.23680.359714.4064×10−54.4154×10−9


### 3.2. Convergence and Error Analysis in the Nonlinear Case

To further verify the accuracy of our proposed method on nonlinear problems, we will calculate a special case of the fourth-order nonlinear system describing fluid flow in this porous channel and compare it with other numerical results.Firstly, to test the convergence of our solution, we define the following residual error function:(35)ResErrF=12j+1∑k=02jF¯Mk2j−F¯M−1k2j2,where F¯M and F¯M−1 denote the mth-order and the m−1th-order wavelet solutions, respectively.Then, we used the numerical solution by means of the software package BVP4c in software Matlab as a reference solutionFref since the original problem did not have an exact analytical solution. Therefore, we define the relative error as follows:(36)ErrM=12j+1∑k=02jF¯Mk2j−Frefk2j2.Based on the definitions of Dauenhauer and Majdalani [12], the axial velocity U and the normal velocity V are given by(37)U=F′R,V=FR.The residual error functions of the axial and normal velocities are therefore expressed as(38)ResErrU=12j+1∑k=02jU¯Mk2j−U¯M−1k2j2,ResErrV=12j+1∑k=02jV¯Mk2j−V¯M−1k2j2.We start our verification forα=0, corresponding to the case of fluid flow through the porous channel with fixed walls. The injection Reynolds number is chosen as R=5, and the convergence-control parameters for U and V are designated as c0=−0.08 and c0=−0.05, respectively. As shown in Figure 3, the residual errors for both U and V reduce very rapidly as the computational order evolves. The least residual error reaches about 10−14 at the computational order M=100 when the resolution level j=6. This figure clearly shows that error reduction can be achieved either by the increment of the resolution level or the increase of the computational order.Figure 3
Residual errors at different resolution levels withc0=−0.08 (for U) and c0=−0.05 (for V). (a) Residual error for U. (b) Residual error for V.
(a)(b)We then check results of the axial velocityU and the normal velocity V for different prescribed points at certain computational orders. As shown in Table 2 and Table 3, the computational results at each discrete point approach to the fixed values very quickly. Furthermore, we compare these results with the numerical ones, and the very excellent agreement is found, as shown in Figure 4. The corresponding computational time and the relative error for different resolution levels are listed in Table 4. It is seen that the computational efficiency is high enough to give highly accurate solutions in one second, and our wavelet solution is very close to the numerical result, and the relative error can reach to 10−11. Note that the wavelet-homotopy computational parameters used here are c0=−0.08 for U and c0=−0.05 for V. The computational order is M=100.Table 2
Different computational order results of the axial velocityU at R=5 and α=0 in the case of j=6 and c0=−0.08.
ηM=1M=10M=20M=50M=80M=10001.5028581.5177731.5231031.5251421.5251681.5251670.1251.4792191.4930601.4979801.4998431.4998651.4998640.2501.4083061.4189361.4226381.4239761.4239851.4239840.3751.2901401.2955521.2972941.2978061.2977981.2977970.5001.1247991.1234761.1227581.1222971.1222771.1222760.6250.9124930.9041200.9009790.8996850.8996660.8996660.7500.6536550.6402570.6356190.6339830.6339780.6339790.8750.3490640.3365400.3325190.3313060.3313140.3313151000000Table 3
Different computational order results of the normal velocityV at R=5 and α=0 in the case of j=6 and c0=−0.05.
ηM=1M=10M=20M=50M=80M=10000000000.1250.1867410.1881340.1889100.1895310.1895880.1895910.2500.3675920.3701730.3716090.3727490.3728520.3728580.3750.5366620.5400300.5418940.5433570.5434840.5434910.5000.6880650.6916460.6936100.6951220.6952460.6952510.6250.8159250.8190470.8207350.8220000.8220940.8220960.7500.9143910.9164440.9175340.9183190.9183700.9183700.8750.9776580.9783920.9787720.9790330.9790470.9790471111111Figure 4
Comparison of the wavelet-homotopy solutions (line) with the numerical results (diamonds and gradients) of (a) the axial velocityU and (b) the normal velocity V for various R at α=0.
(a)(b)Table 4
The computational time and relative error of axial velocityU and normal velocity V at different resolution levels.
Resolution levelj=3j=4j=5j=6UCPU time (sec.)0.22820.27350.31320.3908ErrM2.2991×10−53.3010×10−74.7763×10−92.7685×10−10VCPU time (sec.)0.24070.28470.29390.3857ErrM1.4304×10−51.9613×10−81.7839×10−101.1653×10−11


## 3.1. Validation in the PDE Case

In this section, we will present a 2D linear problem to verify the accuracy and efficiency of our proposed method. The governing equation considered here is as follows:(24)∂2f∂x2+∂2f∂y2+2π2f1=0,where(25)f1x,y=cosπxcosπy.The corresponding homogeneous Neumann boundary conditions are(26)∂f∂xx=0=0,∂f∂xx=1=0,∂f∂yy=0=0,∂f∂yy=1=0.Note that equation (24) has the analytical solution(27)fex,y=cosπxcosπy.In our approach, by modifying the coefficient matrices, the Coiflet wavelet basis is expressed by(28)hj,kx=φj,kxp0,1,i⟶0,p1,1,i⟶0,hj,ly=φj,lyp0,1,i⟶0,p1,1,i⟶0.Furthermore,fx,y can be expressed by Coiflet wavelet as(29)fx,y≈∑k=02j∑l=02jfk2j,l2jhj,kxhj,ly,f1x,y≈∑k=02j∑l=02jfk2j,l2jφj,kxφj,ly.Therefore, equation (24) is expressed as(30)∑k=02j∑l=02jfk2j,l2jd2hj,kxdx2hj,ly+∑k=02j∑l=02jfk2j,l2jhj,kxd2hj,lydy2+2π2∑k=02j∑l=02jf1k2j,l2jφj,kxφj,ly=0.As described above, multiplying both sides of equation (30) by hj,nxhj,myn,m=0∼2j and then applying the multiple integrals on the domain 0,12, we obtain the wavelet-Galerkin equation as(31)A1˜T+A2˜T⋅vecf+B˜T⋅vecf1=0,where vec(⋅) is a straighten vector operator which pulls a matrix into a vector row by row and f,f1 is a 2j+1×2j+1 matrix defined by(32)f=fk,l=fk2j,l2j,f1=fk,l1=f1k2j,l2j,where the coefficient matrices are obtained through(33)A1˜=Γk,nj,2⊗Γl,mj,0,A2˜=Γk,nj,0⊗Γl,mj,2,B˜=2π2⋅Γ¯k,nj,0⊗Γ¯l,mj,0.Here,⊗ denotes the operator of the Kronecker tensor product, A1˜,A2˜,andB˜ are 2j+1×2j+1 matrices, and the generalized connection coefficients Γ and Γ¯ are defined as equation (23).To illustrate the error level of our solution, we define the average square errorErrSQ based on the exact solution (27) as follows:(34)ErrSQ=12j+12∑k=02j∑l=02jfk2j,l2j−fek2j,l2j2.The solution of equation (24) is shown in Figure 2, whose residual error is listed in Table 1. It is found that the average square error decreases rapidly with the resolution level j increasing. For j=7, especially, the error reaches to 10−9 within 1 minute CPU time. This completely demonstrates the accuracy and great efficiency of our proposed method.Figure 2
The wavelet solution of our example problem withj=6.Table 1
The computational time and average square error ErrSQ for equation (24).
Resolution levelj=4j=5j=6j=7CPU time (sec.)0.43490.63841.715332.3233ErrSQ2754.23680.359714.4064×10−54.4154×10−9


## 3.2. Convergence and Error Analysis in the Nonlinear Case

To further verify the accuracy of our proposed method on nonlinear problems, we will calculate a special case of the fourth-order nonlinear system describing fluid flow in this porous channel and compare it with other numerical results.Firstly, to test the convergence of our solution, we define the following residual error function:(35)ResErrF=12j+1∑k=02jF¯Mk2j−F¯M−1k2j2,where F¯M and F¯M−1 denote the mth-order and the m−1th-order wavelet solutions, respectively.Then, we used the numerical solution by means of the software package BVP4c in software Matlab as a reference solutionFref since the original problem did not have an exact analytical solution. Therefore, we define the relative error as follows:(36)ErrM=12j+1∑k=02jF¯Mk2j−Frefk2j2.Based on the definitions of Dauenhauer and Majdalani [12], the axial velocity U and the normal velocity V are given by(37)U=F′R,V=FR.The residual error functions of the axial and normal velocities are therefore expressed as(38)ResErrU=12j+1∑k=02jU¯Mk2j−U¯M−1k2j2,ResErrV=12j+1∑k=02jV¯Mk2j−V¯M−1k2j2.We start our verification forα=0, corresponding to the case of fluid flow through the porous channel with fixed walls. The injection Reynolds number is chosen as R=5, and the convergence-control parameters for U and V are designated as c0=−0.08 and c0=−0.05, respectively. As shown in Figure 3, the residual errors for both U and V reduce very rapidly as the computational order evolves. The least residual error reaches about 10−14 at the computational order M=100 when the resolution level j=6. This figure clearly shows that error reduction can be achieved either by the increment of the resolution level or the increase of the computational order.Figure 3
Residual errors at different resolution levels withc0=−0.08 (for U) and c0=−0.05 (for V). (a) Residual error for U. (b) Residual error for V.
(a)(b)We then check results of the axial velocityU and the normal velocity V for different prescribed points at certain computational orders. As shown in Table 2 and Table 3, the computational results at each discrete point approach to the fixed values very quickly. Furthermore, we compare these results with the numerical ones, and the very excellent agreement is found, as shown in Figure 4. The corresponding computational time and the relative error for different resolution levels are listed in Table 4. It is seen that the computational efficiency is high enough to give highly accurate solutions in one second, and our wavelet solution is very close to the numerical result, and the relative error can reach to 10−11. Note that the wavelet-homotopy computational parameters used here are c0=−0.08 for U and c0=−0.05 for V. The computational order is M=100.Table 2
Different computational order results of the axial velocityU at R=5 and α=0 in the case of j=6 and c0=−0.08.
ηM=1M=10M=20M=50M=80M=10001.5028581.5177731.5231031.5251421.5251681.5251670.1251.4792191.4930601.4979801.4998431.4998651.4998640.2501.4083061.4189361.4226381.4239761.4239851.4239840.3751.2901401.2955521.2972941.2978061.2977981.2977970.5001.1247991.1234761.1227581.1222971.1222771.1222760.6250.9124930.9041200.9009790.8996850.8996660.8996660.7500.6536550.6402570.6356190.6339830.6339780.6339790.8750.3490640.3365400.3325190.3313060.3313140.3313151000000Table 3
Different computational order results of the normal velocityV at R=5 and α=0 in the case of j=6 and c0=−0.05.
ηM=1M=10M=20M=50M=80M=10000000000.1250.1867410.1881340.1889100.1895310.1895880.1895910.2500.3675920.3701730.3716090.3727490.3728520.3728580.3750.5366620.5400300.5418940.5433570.5434840.5434910.5000.6880650.6916460.6936100.6951220.6952460.6952510.6250.8159250.8190470.8207350.8220000.8220940.8220960.7500.9143910.9164440.9175340.9183190.9183700.9183700.8750.9776580.9783920.9787720.9790330.9790470.9790471111111Figure 4
Comparison of the wavelet-homotopy solutions (line) with the numerical results (diamonds and gradients) of (a) the axial velocityU and (b) the normal velocity V for various R at α=0.
(a)(b)Table 4
The computational time and relative error of axial velocityU and normal velocity V at different resolution levels.
Resolution levelj=3j=4j=5j=6UCPU time (sec.)0.22820.27350.31320.3908ErrM2.2991×10−53.3010×10−74.7763×10−92.7685×10−10VCPU time (sec.)0.24070.28470.29390.3857ErrM1.4304×10−51.9613×10−81.7839×10−101.1653×10−11


## 4. Results and Discussion

As discussed above, investigations of the special case of the fourth-order nonlinear system show excellent computational efficiency of our proposed method. Therefore, we apply it to the common cases of the fourth-order nonlinear system for a range of wall expansion ratiosα with different injection or suction Reynolds number R. Without loss of generality, the resolution level j=6 is used for all the following calculations, but the convergence-control parameter needs to be reasonably selected based on the nonlinear strength of the problem at different R and α.It has been known that the traditional homotopy analysis method is very powerful to give highly accurate solutions of nonlinear problems. Here, we compare the computational capability between the traditional HAM technique and our proposed method, firstly. As moderate physical parameters are prescribed, both approaches give very excellent results. However, our proposed technique is superior to the traditional HAM technique for large physical parameters such asα=15 and α=−20, as shown in Figures 5 and 6. In these special cases, the homotopy Padé technique [13] has to be used to optimize the HAM results for high-precision solutions.Figure 5
Comparison of the wavelet-homotopy solutions (lines) with the traditional HAM results (circles) over a range of wall expansion ratiosα at injection Reynolds number R=10.Figure 6
Comparison of our wavelet-homotopy solutions (lines) with the HAM results (circles) over a range of wall expansion ratiosα at suction Reynolds number R=−1.Tables5 and 6 give the comparison of CPU times for the maximum values of the axial velocity U0 by the wavelet-homotopy method and HAM technique with different α for R=10 and R=−1, respectively. Note that the residual errors used here are just an indication to illustrate solution validity. It is seen that our proposed technique has superior computational efficiency over the traditional HAM technique; especially, nonlinearity becomes very strong. For example, in the case of R=10 and α=5, the HAM technique requires 28,450.3 seconds to complete the calculation, but the wavelet-homotopy technique only takes 0.6066 second. Note that the resolution level j=6 is employed in our computation with the convergent-control parameters being listed in Tables 7 and 8.Table 5
Comparison of CPU times forU(0) given by the wavelet-homotopy method and HAM technique with different α in the case of R=10.
αU(0)Residual errorCPU time (s)CurrentHAMCurrentHAMCurrentHAM202.81451—2.2777×10−4—0.7590—152.46663—1.9072×10−8—0.7298—102.12775—1.0711×10−9—0.6611—51.811531.811619.9394×10−114.5521×10−50.606628,450.341.752611.752644.8721×10−111.3808×10−60.530014,106.431.695571.695582.2598×10−112.9759×10−60.50023,506.9221.640631.640641.0360×10−112.3063×10−60.5194980.06611.588021.588035.3207×10−122.4354×10−60.4536292.55001.537961.537963.4777×10−121.3808×10−60.373363.9997−11.490591.490622.5996×10−122.2670×10−70.502932.0802−2.51.425031.425071.3004×10−121.1737×10−60.504632.8491−51.331181.331158.6448×10−123.6382×10−60.6230480.474−101.19967—4.7928×10−10—0.6400—−201.08636—6.8731×10−8—0.7285—Table 6
Comparison of CPU times forU(0) given by the wavelet-homotopy method and HAM technique with different α in the case of R=−1.
αU(0)Residual errorCPU time (s)CurrentHAMCurrentHAMCurrentHAM52.68605—5.5546×10−10—0.7220—2.51.88572—3.0587×10−12—0.6875—21.783801.783973.1057×10−128.7826×10−60.6254974.45011.618021.618094.8021×10−131.5901×10−60.551837.752901.492321.492282.7859×10−154.9944×10−60.56307.22693−11.396861.396632.1015×10−138.4533×10−60.6512986.366−2.51.293811.293844.7729×10−128.9765×10−60.59601,682.30−51.190851.190931.7206×10−101.8075×10−60.66222,181.17−101.10314—4.2413×10−9—0.5935—−201.05330—7.4941×10−9—0.7536—Table 7
The values ofc0 used to calculate U for different R and α.
UR=10R=50R=200αc0αc0αc020−0.0320−0.05100−0.00515−0.05510−0.0650−0.0210−0.075−0.0620−0.030∼5−0.080−0.060−0.05−1∼−5−0.08−5−0.06−20−0.05−10−0.08−20−0.06−50−0.06−20−0.09−50−0.06−100−0.02R=−1R=−5R=−10αc0αc0αc05−0.162.5−0.51−0.162.5−0.142−0.20−0.082−0.131−0.08−1−0.081−0.110−0.08−5−0.080, −1−0.08−1−0.08−20−0.07−2.5−0.02−5−0.08−5−0.013−20−0.08−10, −20−0.01Table 8
The values ofc0 used to calculate V for different R and α.
VR=5R=20R=−1R=−5αc0αc0αc0αc015−0.12520−0.055−0.082.5−0.55−0.055−0.062−0.081−0.180−0.050−0.071−0.150−0.12−5−0.05−5−0.080−0.1−5−0.1−20−0.2−20−0.11−5−0.08−20−0.08−20−0.06Based on the above analysis, we further obtain the solutions for various physical parameters readily by means of our proposed technique, as shown in Figures7 and 8. It can be seen that the values of the axial velocity U reach the highest at the middle section of the porous channel which is η=0, and as the wall expansion ratios, α, increase, the maximum value of the axial velocity becomes gradually larger; however, this situation is reversed near the walls. In addition, increasing wall expansion causes more rapid change in the value of the axial velocity from the walls to the middle section of the porous channel, while increasing wall contraction results in a smoother axial velocity curve with a wider range of plateaus. What is noteworthy is that flow reversal happens near the walls at R=−5, α=2 and 2.5, and R=−10, α=1 cases, which is the result of the rapid expansion of the wall creating a suction area near the wall, so the fluid is forced to move in the direction of the wall to occupy the space created by the wall expansion process. This is also consistent with the findings by Dauenhauer and Majdalani [12].Figure 7
The values of axial velocityU obtained by the wavelet-homotopy method in different wall expansion ratios α at (a) R=50 and (b) R=200.
(a)(b)Figure 8
The values of axial velocityU obtained by the wavelet-homotopy method in different wall expansion ratios α at (a) R=−5 and (b) R=−10.
(a)(b)


## 5. Conclusions

In this paper, the wavelet-homotopy method has been developed to obtain highly accurate solutions of fourth-order nonlinear boundary value problems describing fluid flows in the porous channel with moving walls. By far, now, homotopy-wavelet technique has been investigated for this kind of problems. The full solution procedure has been given in detail. We have performed wavelet approximation on the homogeneous Newman boundary condition by modifying the Coiflet wavelet basis function. The efficiency and accuracy of the Coiflet wavelet method have been verified by an example of a linear two-dimensional problem. Furthermore, the application of our homotopy-based Coiflet wavelet method to nonlinear problems has been carried out. By comparing with the results of the traditional homotopy analysis method, the reliability and efficiency of our Coiflet wavelet-homotopy method in solving strong nonlinear problems are further explained. The solutions of classic problems with liquid flow in porous channels for various physical parameters have been obtained readily by our proposed method. It is expected that this technique is expected to give accurate solutions for more complicated problems with strong nonlinearity.



---

*Source: 5739648-2020-05-15.xml*