Complexity

Volume 2018, Article ID 8269541, 11 pages

https://doi.org/10.1155/2018/8269541

## A Numerical Algorithm for Solving Higher-Order Nonlinear BVPs with an Application on Fluid Flow over a Shrinking Permeable Infinite Long Cylinder

^{1}Department of Physics, United Arab Emirates University, P.O. Box 15551, AlAin, UAE^{2}Department of Mathematical Sciences, United Arab Emirates University, P.O. Box 15551, AlAin, UAE

Correspondence should be addressed to Qasem M. Al-Mdallal; ea.ca.ueau@lalladmla.q

Received 5 November 2017; Accepted 12 February 2018; Published 12 March 2018

Academic Editor: Mohamed Belhaq

Copyright © 2018 Laila Y. Al Sakkaf et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. Introduction

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 [34]. 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 [49].

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. [8] and Elnajjar et al. [10]. 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 [10] 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 [10] 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

#### 5. Validation

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.