A Numerical Algorithm for Solving Higher-Order Nonlinear BVPs with an Application on Fluid Flow over a Shrinking Permeable Infinite Long Cylinder
We present an efficient iterative power series method for nonlinear boundary-value problems that treats the typical divergence problem and increases arbitrarily the radius of convergence. This method is based on expanding the solution around an iterative initial point. We employ this method to study the unsteady, viscous, and incompressible laminar flow and heat transfer over a shrinking permeable cylinder. More precisely, we solve the unsteady nonlinear Navier–Stokes and energy equations after reducing them to a system of nonlinear boundary-value problems of ordinary differential equations. The present method successfully captures dual solutions for both the flow and heat transfer fields and a unique solution at a specific critical unsteadiness parameter. Comparisons with previous numerical methods and an exact solution verify the validity, accuracy, and efficiency of the present method.
Numerous phenomena in engineering and applied science fields are governed by nonlinear boundary-value problems (BVPs). Therefore, BVPs have received a huge attention from mathematicians, physicists, and engineers for the sake of finding and analyzing their solutions. Generally speaking, finding the analytical solutions for nonlinear BVPs is far from trivial and often is impossible. Therefore, many numerical techniques have been developed to solve such type of problems. These methods include Adomian’s decomposition method, homotopy perturbation method, variational iteration method, optimal homotopy asymptotic method, operational matrices techniques based on various orthogonal polynomials and wavelets, finite difference method, and spectral methods; the reader is referred to [1–6] and references therein.
The fluid dynamics and heat transfer of a viscous incompressible fluid flowing past stretching surfaces, such as a sheet or tube, have attracted considerable interests of many researchers because of their importance in many industrial applications such as the quality of certain products. One of the most interesting conditions for stretching surfaces problems is the velocity at the surface, where it mainly figures the characteristics of the fluid based upon two essential factors, fluid viscosity and suction parameter. A remarkable interest of several researchers concentrated on tracking the existence of dual solutions for the flow within a certain range of unsteadiness and suction parameters [7–17]. Although, the literature reveals numerous research papers discussing the flow over a stretching sheet and moving plate [18–26], there are only few studies focusing on the problem of flowing past a stretching cylinder or tube; see [8–10] and references therein.
It is established that the differential equations describing the fluid flow BVPs are highly nonlinear and demand extremely accurate numerical schemes. In this work, we show that an iterative procedure, based on successive power series expansions, provides one such high accuracy numerical scheme. Often, using power series solutions turns to be useless because the resulting solution diverges at a finite radius of convergence. The divergence is intrinsic to the nature of the solution in the sense that it persists to exist even with an infinite power series expansion. The method presented here solves this problem showing that the radius of convergence can be delayed arbitrarily to any large value. This value could, in principle, approach infinity achieving exact solutions.
Among the many different methods of solving nonlinear differential equations [27–33], the power series method is the most straight forward and efficient . It has been used as a powerful numerical scheme for many problems [35–43] including chaotic systems [44–47]. Many numerical algorithms and codes have been developed based on this method [34–36, 44–48]. However, the abovementioned finiteness of radius of convergence is a serious problem that hinders the use of this method to wide class of differential equations, in particular the nonlinear ones. For instance, the nonlinear Schrödinger equation (NLSE) with cubic nonlinearity has as a solution. Using the power series method to solve this equation produces the power series of , which is valid only for .
Recognizing that a powerful numerical scheme based on this method is already established [34–36, 44–48], we nonetheless present a thorough investigation of the error associated with this method with the aim of showing how we can systemically reduce errors to infinitesimal values while having the Central Processing Unit (CPU) time within a reasonable range. We will show robustness and efficiency of the method via a highly demanding fluid flow boundary-value problem. Therefore, solving the problem of finite radius of convergence will open the door wide for applying the power series method to much larger class of differential equations, particularly the nonlinear ones.
Briefly, the present technique is based on iterative power series expansions of the solution. The domain of the independent variable, say , is divided into a number of segments each of width , where is smaller than the radius of convergence. A power series solution is obtained by expanding the solution around the left end of the first segment using the initial conditions given with the problem. Similarly, a power series solution is obtained by expanding around the start of the second segment but now using the first series to calculate the initial conditions. This is repeated times till a solution at is obtained. In the limit and the series solution becomes an exact solution. This scheme is effectively equivalent to an iterative procedure of repeated iterative calculation of the recursion relations of the power series in the first segment. Another aim of this paper is to apply this method to study the unsteady flow and heat transfer characteristics of fluid flow over a shrinking permeable infinite long cylinder. We will show that the iterative numerical scheme resulting from this method is exceeding the efficiency of typical numerical methods used. In addition, we managed to find an exact solution which enabled us to calculate accurately the error for a finite value of number of iterations . It should be noted that the present work is a part of the Master thesis .
The rest of the paper is organized as follows. In Section 2, a mathematical representation of our method is illustrated using a general form of nonlinear ordinary differential equation, while henceforth we call it iterative power series method. Sections 3 and 4 display the implementation of the iterative power series method on the heat and mass transfer model. Section 5 focuses on analyzing the validity of this present technique and demonstrating its efficiency by drawing comparisons with the achieved exact analytical solution and other numerical methods. In Section 6, we analyze and discuss the properties of the solutions obtained. We end with a summary of our main conclusions in Section 7.
2. General Scheme of the Iterative Power Series (IPS) Method
In this section, we give a brief description of the present method. Consider a general ordinary differential equation of the formwith initial conditionswhere is the th derivative of , are real constants, and is a known function. The factor is introduced, without loss of generality, for the constants to correspond to the coefficients of the power series expansion below. At first, we divide the interval into a number of identical segments each of width . Then we expand in a power series around the beginning of each interval, namely,where is the power series expansion around the start of the th segment and are the coefficients of the power series. Recursion relations between the coefficients are obtained upon substituting the power series solution, (3), into the differential equation, (1), which can be expressed in terms of the first coefficients corresponding to the initial conditionswhere denotes the set of coefficients . The essential idea of the IPS method is to calculate the coefficients of the th power series from the th series according to (2),which upon using (3) readsand is simplified toand then imposes the condition . Here, is the binomial function. The last equation is the basis for the IPS algorithm. Starting from the initial conditions for the power series of the zeroth interval, an iterative application of (7) leads to the coefficients of the th interval, namely, which give the solution at the desired point, ,Both analytical and numerical schemes may be deduced from this algorithm. For the numerical scheme, the value of used is inserted as a number . On the other hand, leaving as a variable results in an analytical solution in terms of a power series in which is equivalent to a functional transformation on the zeroth order series; that is, the coefficients of the th series are functional transformation of the th series. In such a case the last power series for the th interval corresponds to such that functional transformations and all power series expansions of the zeroth up to th intervals will be included in the th expansion.
The coefficient of each th expansion represents the value of the solution at , which gives a discrete representation of . Therefore, in the limit , the discrete representation turns to a continuous one and thus we conjecture that the exact solution is obtained in the limits of and
3. Heat and Mass Transfer Model
In this section we will employ the present numerical technique on heat and mass transfer over a shrinking permeable cylinder described in Zaimi et al.  and Elnajjar et al. . For completeness, we redescribe precisely the physical model. The flow is considered an unsteady, laminar, viscous, and incompressible fluid with uniform velocity and uniform temperature over a permeable shrinking circular cylinder. The cylinder is assumed to be infinitely long and the flow has constant properties. The diameter of the cylinder is assumed to be time dependent with the radius , where is a positive constant, is the constant of expansion/contraction strength, and is the time. Clearly, the cylinder’s radius is shrinking with time if is positive and stretching with time if is negative. Notice that since the flow is axisymmetric, the flow field should be a function of the radial coordinate, , and the longitudinal. The governing equations for the unsteady and incompressible fluid without body force are the continuity, momentum, and energy equations. These equations in cylindrical coordinate system, , are given bywhere and are the polar coordinates in the radial and axial directions, respectively, and are the fluid velocity components in the radial and axial directions, respectively, and is the fluid temperature. The function represents the fluid pressure and the parameters , , and denote the fluid viscosity, the fluid density, and the fluid thermal diffusivity, respectively. Notice that we assumed that there is no azimuthal velocity component. The assumed boundary conditions associated with (10) for the velocity components and the temperature are given bywhere is the constant surface temperature and is a positive constant.
The similarity transformations which convert (10) into nonlinear ordinary differential equations are given by [8, 10]where and is the similarity variable given byIn addition, it should be noted that represents the dimensionless stream function and represents the normalized temperature. Applying the above similarity transformations, (10) and the boundary conditions (11) reduce tosubject towhere is the unsteadiness parameter representing the strength of contraction/expansion, is the suction parameter, and is the Prandtl number.
Our main target is to solve (14) and (15), subject to the boundary conditions (16), using the present technique in the ranges and at . In this study, we will analyze the normalized skin friction coefficient, , and the normalized heat transfer rate, .
4. IPS Method for Heat and Mass Transfer Model
The following is a detailed implementation of the IPS method used to solve (14) and (15) subject to the boundary conditions (16). It is worth mentioning that (14) includes only the variable , while (15) includes both and . Therefore, it is more convenient to find the variable from (14) and then solve (15). Furthermore, for the sake of simplicity, we render (14) to an initial-value problem; that is,withwhere must be chosen using the shooting method  so that the solution satisfies the boundary condition . Notice that by setting different initial values for in the shooting method, dual solutions are obtained.
As mentioned in Section 2, we start by expanding in power series around the initial point and the derivatives, , , and can be calculated simply by differentiating this series. Substituted in (14), the coefficients, , for 3, can be found recursively in terms of the initial conditions , , and through the recursion relations. The first two recursion relations are given byRecalculating , , and at givesNow, , , and play the role of the initial conditions for the next series expansion, where we expand the solution and its derivatives in power series around Resubstituting these power series expansions in the differential equation, we get the new recursion relations . The next iterative step is to calculate and its first two derivatives at which will give the initial conditions for the new power series. Repeating this process times, the general forms of the first two recursion relations of (14) are found to bewhereIn summary, the IPS procedure can be reduced to the following algorithm:where and are the recursion relations obtained from the differential equation. We have removed the superscripts that indicate the index of the iteration for convenience. The scheme is thus described simply as follows: one starts with (26) to calculate , followed by updating the initial conditions according to (27) and then using the updated values back in (26) and so on. The procedure has to be repeated times with .
Similarly, for the energy equationsubject towhere must be chosen using the shooting method  so that the solution satisfies the boundary condition . We expand in power series around the same initial pointand similarly for the derivatives, and . Substituted in (15), the coefficients, , for 2, are found recursively in terms of the initial conditions and through the recursion relations. Employing the IPS method, the first two recursion relations are found to be where
In this section we aim at demonstrating the performance and efficiency of the present numerical scheme. Firstly, in order to obtain accurate numerical results, we have to pay attention to the selection of the numerical algorithm parameters, , , and . To achieve this target we focus on studying the following:withwhere is updated using the shooting method to satisfy the condition . Notice that several recent studies such as [10–16] reported the existence of a critical value of (named ) at which the problem has no solution (for ), only one solution (at ), and dual solutions (for ). It is worth mentioning that we succeeded in finding the explicit analytical form of the first solution for (33) subject to (34) under a condition ; that isThis exact solution will play a crucial role in proving the advantages of our numerical scheme.
The approximate solutions of the problems (33) and (34) at three different iterations , and , together with the exact solution obtained by (35) when and , are displayed in Figure 1. It is clearly seen that increasing the number of iterations in the IPS method delays the divergence point.
To achieve an “optimal choice” of , we solve the problem with . Table 1 shows the values of up to 50 digits corresponding to the values of . It is clearly seen that the value of stabilizes at around ; hence we choose as the optimal value for the rest of the calculations in the entire paper. It should be noted that most of the used numerical schemes for such type of problems, [8–10, 17], had chosen to represent the infinity which, definitely, gave lower order of accuracy. This conclusion can be easily tested via the exact solution (35) which gives and . However, we will only show the interval up to for the rest of the coming figures.
The upper bound of the error in the IPS method can be estimated as follows. At each iterative step an error of order results from terminating the power series at . This error will be magnified times due to the iterative procedure. As a result, the upper bound of the error of the IPS method is
Table 2 presents the CPU time in seconds, which is machine-dependent, versus the upper bound of the error, , of the IPS method for the first solution at and , where varies from 3 to 8.
The exact solution, (35), provides a unique possibility of calculating the error of the IPS method and comparing it with that of other numerical methods. Figure 2 presents a comparison between the error of the IPS method and the explicit Runge–Kutta method of order four (ERK4) for the problem at and . The advantages of the present technique over the other one are notable.
A further comparison is done in Figure 3 on the normalized skin friction coefficient, , as a function of with Zaimi et al.  and Elnajjar et al.  for the case when . Excellent agreements are obtained. It should be mentioned herein that Elnajjar et al.  used a combination of the implicit Runge–Kutta method and the shooting method while Zaimi et al.  implemented the shooting method described in the book by Jaluria and Torrance .
It is worth mentioning here that our technique successfully showed its definite capability to exceed the machine precision which is . A clear sign on this is the systematic reduction in the error in Table 2 when compared with the exact solution. Moreover, all computations are performed with at least 14 decimal digits of precision, knowing that all computations are operated using Mathematica software 10.4 and carried out on a Lenovo PC with the following specifications: model: Z470, processor: Core(TM) i5-2430M CPU @ 2.40 GHz, system type: 64-bit, and installed memory (RAM): 4.00 GB.
6. Results and Discussion
In this section, we discuss the effects of both suction parameter, , and unsteadiness parameter, , on the velocity profile, , the normalized skin friction coefficient, , the temperature profile, , and the heat transfer rate, . The numerical simulations are conducted at a fixed Prandtl number, , while the ranges considered for the other parameters are and .
Figure 4 shows the first and second solutions of the velocity profiles for with a fixed unsteadiness parameter, . It is clearly seen that the first solution for the fluid velocity inside the boundary layer region increases as increases, while the second solution shows an opposite trend. In addition, the two solutions of the velocity profile become steeper with higher magnitudes as increases. These observations emphasize the effect of increasing the suction parameter of the cylinder’s wall which is to decrease the boundary layer thickness. Consequently, increasing the suction parameter causes an increment in the normalized skin friction coefficient for the first solution and decrement in the normalized skin friction coefficient for the second solution, as clearly shown in Figure 5. These findings are consistent with the results reported by Elnajjar et al.  and Zaimi et al. .
Figure 6 presents the temperature profiles of the fluid flow, , at and . It is obviously noticeable that both solutions for temperature profiles admit similar behaviour, where they become wider and more relaxed as the suction parameter decreases. These kinds of behaviour inspire us to conclude that the developed thermal boundary layer and the corresponding rate of heat transfer are decreasing as increases. However, the second solution depicts more relaxed behaviour compared with the first solution. This slight difference between the first and second temperature profiles indicates that the second solution reflects higher thermal boundary layer than the first solution and, thus, a larger rate of the heat transfer as confirmed by Figure 7.
The variation of both the normalized skin friction coefficient, , and the heat transfer rate, , as functions of , are shown, respectively, in Figures 5 and 7 for . The results demonstrate the existence of a critical value in the -domain at which the problem has no solution for , only one solution at , and dual solutions for . Figure 5 shows that increases as increases which is due to the increase in the surface shear stress coefficient. Moreover, we observe that is decreasing with . However, Figure 7 clearly shows that increasing will definitely increase the heat transfer rate while increasing causes a decrease in the heat transfer rate.
Figure 8 displays the first and second solutions of the velocity profiles for with a fixed value of the suction parameter, . Generally speaking, the behaviour of is very similar to the case of the variable suction parameter; that is, increasing the unsteadiness parameter produces steeper behaviour in the velocity profiles for the first solution while the second solution shows an opposite trend. In agreement with the case of the variable suction parameter in Figure 4, increasing the unsteadiness parameter will then cause a reduction in the thickness of the boundary layer.
Figure 9 presents the temperature profiles of the fluid flow, , at and . Clearly, the increase in the unsteadiness parameter or the suction parameter leads to the same trend.
Finally, we end up our discussion with Figure 10 which presents an overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the domain for and . The straight line in this figure represents the occurrence of unique solution of the problem.
In this work, we presented a numerical technique for solving nonlinear BVPs based on iterative power series solutions. We have demonstrated its efficiency and accuracy through validation against the numerical ERK4. We have shown that our method excels over the ERK4 by orders of magnitude in accuracy. Moreover, the accuracy in our method is systematically controlled such that the error can be reduced to any arbitrary small value. We successfully studied the unsteady viscous flow over a contracting cylinder using the present technique. The velocity and temperature profiles of the ordinary version of Navier–Stokes equations for different suction and unsteadiness parameters are calculated. We have obtained the exact solution of the first solution of the fluid flow under a specific condition, , and have employed it to emphasize the efficiency of the present numerical technique.
The convergence analysis of the IPS method is considered for a future work. However, we strongly believe that the method is convergent. This is conjectured by the systematic reduction in error upon increasing the number of iterations or the number of terms in the power series. We believe that this technique will serve researchers in different fields working on nonlinear systems. In particular, the technique will be very useful for systems described by nonintegrable nonlinear differential equations.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
The authors would like to express their sincere appreciation to the United Arab Emirates University, Al Ain, UAE, for providing the financial support with the UPAR (7) 2015 and the UPAR (4) 2016 grants.
Y. Jaluria and K. E. Torrance, Comput Heat Transfer, 2003.
A. D. Polyanin and V. F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, CRC press, 2004.View at: MathSciNet
M. A. Salle and V. B. Matveev, Darboux transformations and solitons, 1991.
L. Y. Al Sakkaf, Iterative power series solutions of nonlinear partial differential equations with applications in nonlinear Science [M.S. thesis], United Arab Emirates University, Al Ain, United Arab Emirates, 2018.