Semi-Analytical Solution of Two-Dimensional Viscous Flow through Expanding/Contracting Gaps with Permeable Walls

: In this research, the analytical methods of the differential transform method (DTM), homotopy asymptotic method (HAM), optimal homotopy asymptotic method (OHAM), Adomian decomposition method (ADM), variation iteration method (VIM) and reproducing kernel Hilbert space method (RKHSM), and the numerical method of the ﬁnite difference method (FDM) for (analytical-numerical) simulation of 2D viscous ﬂow along expanding/contracting channels with permeable borders are carried out. The solutions for analytical method are obtained in series form (and the series are convergent), while for the numerical method the solution is obtained taking into account approximation techniques of second-order accuracy. The OHAM and HAM provide an appropriate method for controlling the convergence of the discretization series and adjusting convergence domains, despite having a problem for large sizes of obtained results in series form; for instance, the size of the series solution for the DTM is very small for the same order of accuracy. It is hard to judge which method is the best and all of them have their advantages and disadvantages. For instance, applying the DTM to BVPs is difﬁcult; however, solving BVPs with the HAM, OHAM and VIM is simple and straightforward. The extracted solutions, in comparison with the computational solutions (shooting procedure combined with a Runge–Kutta fourth-order scheme, ﬁnite difference method), demonstrate remarkable accuracy. Finally, CPU time, average error and residual error for different cases are presented in tables and ﬁgures.


Introduction
Significant interest has been paid to the analysis of nonlinear challenges in different fields of nature and engineering. To solve this type of problem, the wide group of analytical solutions and computational simulations has been considered. A general analytical method to solve the non-linear differential equations is the perturbation technique [1,2]. Practical

Differential Transform Method
Using [9,12,13], the major description of the differential transformation can be described easily. Consider a function u(x) that is defined by a power series with a center in a point x 0 . The differential reduction of u(x) can be defined as Here U(k) is the transformed function. In addition, the inverse transformation is determined as In combination with Equations (1) and (2), the following equation is generated: Considering Equation (3), Taylor series expansion is the basis of differential transformation. This technique does not allow us to perform an evaluation of the derivatives symbolically. In physical problems, the finite series of u(x) is considered and Equation (2) can be formulated as assuming that u(x) = ∑ ∞ k=m+1 (x − x 0 ) k U(k) has negligible values. Mainly, the magnitude of m is defined by convergences of the series parameters.

Solution
In this section, six analytical methods and two numerical techniques are applied to work out the boundary value challenge. Hence, 2D viscous circulation in a channel bordered by two moved porous plates can be described using the following equations ( Figure 12 shows streamlines and velocity field for Re = 5 and α = 1): with additional border restrictions where α is a dimensionless border stretching rate and Re is the penetration Reynolds number.

Differential Transform Method
Employing the differential conversion from Equation (17), the following form of equation is obtained: where F(k) denotes the differential transforms of n f (y) and it is defined as The reduced border restrictions are where a and b are constants. Employing Equation (22) for Equation (20) and by the recurrent technique, all magnitudes of F(k) are calculated. Hence, substituting all F(k) into Equation (21), the series of solutions is defined as below: Then, by applying the boundary conditions of Equation (19), the constants a and b are calculated. For example, for Re = 1 and α = 1, the values and the solutions are as follows: +1.607412616389382y − 0.732958779309001y 3 +0.1465917558618002y 5 − 0.024487625417451822y 7 +0.0040548656166754194y 9 − 0.0007243900882554326y 11 +0.00014026681559397154y 13 − 0.0000287098687570248y 15 (25) For Re = 1 and α = −1 the values and the solutions are as follows:
(1 − p) ∂ 4 ϕ(y,p) The f (y) and H(p) are considered as follows: By using f (y) and H(p) from Equation (32) for Equation (30) and some transformation with the help of powers of p-terms, the following solutions can be obtained.
p 0 : Zero-order problem from which we obtain f 0 (y) = 1 2 3y − y 3 (34) it is obtained that and therefore the terms f 3 (y) are too large to be illustrated graphically. As a result, the final form of f (y) is Using Equation (14) for f (y) in Equation (17), R(y,C 1 ,C 2 ,C 3 ) and J are defined as follows: The constants C 1 , C 2 and C 3 are obtained as follows: For example, for the case of Re = 1 and α = 1, the following 3rd-order approximate solution is determined.

Homotopy Analysis Method
Dinarvand and Rashidi [26] studied this problem by the HAM. Yabushita et al. [23] employed the HAM to solve two coupled nonlinear ODEs by an optimization technique. The optimal characteristics are calculated by minimization of the square residual error defined in the considered physical domain. In the present research, the same as for the OHAM (Equations (14)-(16)), the optimal convergence-control parameter ( ) is found below.
The nth-order approximation of the solution defined using the HAM f (y) is which depends on the convergence characteristic ( ). Let denote the square residual error of the control relation. The optimal value of ( ) is gained by solving the nonlinear algebraic equation.

Adomian Decomposition Method
In this section, the ADM is employed for the solution of Equation (17) along with the boundary conditions (18) and (19). Let us introduce the fourth-order derivative operator L and inverse operators L −1 as follows: Thus Equation (17) becomes The function f (η) is The remaining terms of (17) can be expressed as Here, a recursive formula is used to find all the components. The exact solutions of (17) are given by Therefore, the RHSs of Equation (50) are given by From (18), (19) and invoking the border restrictions The solutions of Equation (50) can therefore be written as The unknowns p, q should be evaluated numerically. Utilizing Equation (55), the initial imposed solutions along with higher-order recursive solutions are Therefore, the solution can be expressed as

Variation Iteration Method
The illustration of the basic concept of the approximate analytical solution, i.e., the variation iteration method (VIM) is as follows.
Let us analyze the standard view of the nonlinear differential relation as The linear function is represented as L, the nonlinear part is N and g(t) is an inhomogeneous part. The construction of the correct operator employing the VIM is presented as Here u 0 (t) is the initial assumed solution obtained using the initial conditions and the suitable choice of unknown initial conditions, u n (ξ) is the restricted function and λ(ξ) is the Lagrangian multiplier obtained using the Wronskian conditions.
Employing the aforesaid methodology, the present problem can be expressed as The Lagrangian multiplier is With the initial choice of the initial approximate solutions along with the final iterative solution are presented.

Reproducing Kernel Hilbert Space Method
The reproducing kernel Hilbert space method (RKHSM) without the Gram-Schmidt orthogonalization process, is analyzed in this manuscript. Some applications of the RKHSM to solve various differential equations are presented in [25][26][27][28]. We consider the reproducing kernel space W 5 [0,1] for Equations (17)- (19), such that f (4) (y) is absolutely continuous and f (5) The inner product and norm are as follows: ) is a non-linear func-tion and F(y) is a right-hand-side function. It is not difficult to understand that L is a bordered linear function. Assume L:W 5 [0,1]→W 1 [0,1], and r y (x) is a reproducing kernel operator for W 1 [0,1]. We choose a dense set {x i } ∞ i=1 in [0,1] and define ϕ i (y) = r x (y) x=x i and ψ i (y) = L * ϕ i (y) where L * is an adjoint function of L and L -1 exists. This can be proved by ξ i (y) = R x (y) x=x i , which are a complete function system in W 5 [0,1]. In this method we assume is an exact solution of Equation (17) Consider an approximate solution for Equation (17) where n is the number of collocation points on [0,1] and ξ i (y) = R x (y) x=x i . One can obtain c i by solving the following iterative system of algebraic equations

Finite Difference Method for ODEs
The main idea of the finite difference method is to approximate the derivatives in ODEs using the finite differences, taking into account the Taylor series and after that solving the obtained system of equations. For this purpose, the uniform mesh is introduced as a first step of this technique. As a result, the numerical solution of the boundary value problem (17)- (19) has been performed using the finite difference method and a uniform mesh. The uniform mesh has been introduced as follows: Derivatives presented in Equation (17) have been approximated using the following central differences of second-order accuracy Using these mentioned central differences, Equation (17) has been reduced to the set of equations that was determined by the successive over-relaxation procedure as follows: The convergence condition used is It should be noted that boundary conditions (18) and (19) have been approximated as follows: The mesh sensitivity analysis has been carried out for this technique using two different uniform meshes of 21 nodes and 40 nodes. Tables 1-4 show that a mesh of 21 nodes is enough for this analysis, because the differences between these two grids are too small.

Finite Difference Method for PDE
It should be noted that the formulated boundary value problem (17)- (19) for ODEs was obtained using the original Navier-Stokes equations with appropriate initial and boundary conditions that were written for laminar, isothermal and incompressible fluid circulation between two horizontal permeable plates.
with additional restrictions that were written as These formulated boundary value problems for partial differential equations were reduced using the stream function and vorticity that were determined as Taking into account these new functions, governing Equations (81)-(83) were rewritten as Additional restrictions (84) and (85) were formulated as Taking into account the following dimensionless parameters
Using these parameters, the non-dimensional governing equations are Non-dimensional restrictions are Here α = ρ·a· . a µ and Re = ρ·a·V w µ . These formulated boundary-value problems (92)-(95) were solved by the finite difference method with second-order accuracy for a uniform mesh. The employed schemes of the second-order were used for the diffusive part and convective part, where the first-order approximation was applied for the unsteady term. The convective terms were discretized by means of the monotonic Samarskii procedure and the central difference scheme was used for diffusive terms.
The solution of Equation (93) was obtained using the Samarskii locally one-dimensional technique with an additional time level: where i and j are the indexes for the nodes over the x and y coordinates, n is the time zone index, ∆τ is the time step, and h x and h y are the grid steps for the x and y coordinates. At the final level of solution, the discretized Equations (96) and (97) were solved by the Thomas method. The difference relation for the stream function was worked out employing the successive over-relaxation procedure as follows where k is an index of iterations and κ is the relaxation parameter. It should be noted that this formulated finite difference technique has previously been validated using the mesh sensitivity analysis and some model problems.

Results and Discussion
The program code was developed by MATHEMATICA and the simulations were defined using a PC with 756 MB of RAM and 2.40 GHz CPU. In the case of the FDM for ODEs and PDEs the computational code was developed using C++ programming language. Tables 1-4 show a comparison of the five analytical methods and the FDM and show the absolute error for computational results and the analytical solutions. We have deliberately provided many digits for comparison in the tables to better show the accuracy of each method. In fact, few digits cannot show the accuracy clearly. The numerical solution, the OHAM and HAM by 3rd-order approximate solution and the DTM (m = 10) by 10th-order approximate solution for the case of Re =5 and α = 1 are presented in Table 1. In addition, for the same case we compared the computational results with the OHAM and HAM by 4th-order approximate solution and the DTM (m = 20) by 20th-order approximate solution. Numerical results obtained by the OHAM and HAM (5th-order) and the DTM (m = 30) for the same case are shown in Table 3. For all approaches, the absolute error diminishes quickly with the order of approximation increment. Table 4 compares the numerical simulation results with the ADM and VIM. In the case of the FDM for ODEs, it should be highlighted that the uniform mesh of 21 nodes that was used was enough for obtaining the convergent result.
From Tables 1-4 it can be highlighted that the OHAM approximate relation defined in this research allows us to obtain better outcomes than the DTM and HAM approximations, and the HAM approximate expression is more accurate than the DTM approximation for this problem. In addition, the ADM is more accurate than the VIM.
The function uc/x (f (y)) obtained by the OHAM, HAM and the DTM at various orders of approximation for the case of contraction combined with suction (Re = −5 and α = −1) are compared with the computational outcomes in Figures 1-3.
The comparison of CPU time given by the OHAM, DTM and HAM at various orders of approximation for the case of Re = 5 and α = 1 is given in Table 5 and Figure 4. Since the DTM method uses the transformed function to convert the differential equations to algebraic equations, it is easier to calculate than the HAM and OHAM. In the other two methods, the differential equations are solved in each iteration. From Table 5 and Figure 4, it is inferred that the OHAM has more unknown parameters and takes more CPU time to converge, especially for high-order approximations.
error diminishes quickly with the order of approximation increment. Table 4 compares the numerical simulation results with the ADM and VIM. In the case of the FDM for ODEs, it should be highlighted that the uniform mesh of 21 nodes that was used was enough for obtaining the convergent result.
From Tables 1-4 it can be highlighted that the OHAM approximate relation defined in this research allows us to obtain better outcomes than the DTM and HAM approximations, and the HAM approximate expression is more accurate than the DTM approximation for this problem. In addition, the ADM is more accurate than the VIM.
The function uc/x (f′(y)) obtained by the OHAM, HAM and the DTM at various orders of approximation for the case of contraction combined with suction (Re = −5 and α = −1) are compared with the computational outcomes in Figures 1-3.     The comparison of CPU time given by the OHAM, DTM and HAM at various orders of approximation for the case of Re = 5 and α = 1 is given in Table 5 and Figure 4. Since the DTM method uses the transformed function to convert the differential equations to algebraic equations, it is easier to calculate than the HAM and OHAM. In the other two methods, the differential equations are solved in each iteration. From Table 5 and Figure 4, it is inferred that the OHAM has more unknown parameters and takes more CPU time to converge, especially for high-order approximations.    The average error Equation (73) for K = 100 obtained by different approaches at different orders of approximation for the case of (Re = 5 and α = 1) is shown in Figure 5. Figure 5 illustrates that the accuracy of the solution obtained by the OHAM is better than the HAM and DTM in all orders of approximation (3 rd order OHAM is better that 3 rd order HAM and DTM with 10 terms of polynomials, …). The average error Equation (73) for K = 100 obtained by different approaches at different orders of approximation for the case of (Re = 5 and α = 1) is shown in Figure 5. Figure 5 illustrates that the accuracy of the solution obtained by the OHAM is better than the HAM and DTM in all orders of approximation (3rd order OHAM is better that 3rd order HAM and DTM with 10 terms of polynomials, . . . ). Substituting the approximate solution into Equation (17) yields the residual error. In Figure 6; Figure 7, we present the residual error of different approaches for the case of Re = 5, α = 1. Figure 6 presents the residual error for the OHAM and HAM by 4th-order approximate solution and the DTM by m = 20, and Figure 7 presents the residual error for the OHAM and HAM by 5th-order approximate solution and the DTM by 30th-order approximate solution (with a polynomial by 30 terms). It can be seen from Figure 6; Figure 7 that the accuracy of the solution defined using the OHAM is very good in both cases. In addition, the residual error obtained by the DTM increases as y increases. This is because in the DTM the Taylor expansion is used about (y = 0). For the case of the ADM and VIM solutions, Figure 8; Figure 9 show the f(y) solution with various Re and α. The ADM shows different results in larger y and especially in high Reynolds numbers. Substituting the approximate solution into Equation (17) yields the residual error. In Figures 6 and 7, we present the residual error of different approaches for the case of Re = 5, α = 1. Figure 6 presents the residual error for the OHAM and HAM by 4th-order approximate solution and the DTM by m = 20, and Figure 7 presents the residual error for the OHAM and HAM by 5th-order approximate solution and the DTM by 30th-order approximate solution (with a polynomial by 30 terms). It can be seen from Figures 6 and 7 that the accuracy of the solution defined using the OHAM is very good in both cases. In addition, the residual error obtained by the DTM increases as y increases. This is because in the DTM the Taylor expansion is used about (y = 0). For the case of the ADM and VIM solutions, Figures 8 and 9 show the f (y) solution with various Re and α. The ADM shows different results in larger y and especially in high Reynolds numbers.

Error
In the RKHSM, the residual error (RError) is calculated using in Table 6, and we consider M = 10 and n = 40 for Figures 10 and 11. The obtained results show the high accuracy of this method for large Reynolds numbers. Figure 12 shows the streamlines and velocity field for Re = 5 and α = 1 that were obtained employing the FDM for PDEs.

Conclusions
In this study, different analytical/numerical methods were proposed to solve 2D viscous fluid motion between expanding/contracting horizontal permeable plates. The results were compared with the HAM, ADM, VIM, … and numerical results (shooting procedure combined with a Runge-Kutta fourth-order integration technique and the FDM), and obtained satisfactory outcomes. The result can be further improved by increasing the order of approximation. These approaches are simple to use because, unlike other numerical and approximate methods, they do not require linearization, discretization or perturbation. The proposed procedures are also valid for nonlinear differential equations. The OHAM and HAM provide a convenient method to manage the convergence and the user can easily adjust the desired convergence regions. It is clear that the OHAM has a high accuracy compared to the DTM. However, having many unknown constants is time-consuming when simulating the square residual errors for high-order approximation, as shown in Figure 4 and Table 5. The HAM uses less CPU time compared to the OHAM, which first has the unique convergence managing characteristic ( )  as an unknown and then gives its effective magnitude by minimizing the square residual error.

Conclusions
In this study, different analytical/numerical methods were proposed to solve 2D viscous fluid motion between expanding/contracting horizontal permeable plates. The results were compared with the HAM, ADM, VIM, . . . and numerical results (shooting procedure combined with a Runge-Kutta fourth-order integration technique and the FDM), and obtained satisfactory outcomes. The result can be further improved by increasing the order of approximation. These approaches are simple to use because, unlike other numerical and approximate methods, they do not require linearization, discretization or perturbation. The proposed procedures are also valid for nonlinear differential equations. The OHAM and HAM provide a convenient method to manage the convergence and the user can easily adjust the desired convergence regions. It is clear that the OHAM has a high accuracy compared to the DTM. However, having many unknown constants is time-consuming when simulating the square residual errors for high-order approximation, as shown in Figure 4 and Table 5. The HAM uses less CPU time compared to the OHAM, which first has the unique convergence managing characteristic ( ) as an unknown and then gives its effective magnitude by minimizing the square residual error.