Journal of Applied Mathematics
Volume 2009 (2009), Article ID 524307, 17 pages
doi:10.1155/2009/524307
Research Article

An Application of Homotopy Analysis to the Viscous Flow Past a Circular Cylinder

Department of Mathematics, University of Benin, Benin-City, P.M.B 1154, Nigeria

Received 19 July 2008; Accepted 6 April 2009

Academic Editor: Bernard Geurts

Copyright © 2009 E. O. Ifidon. 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 consider the application of a new analytic method based on homotopy analysis to the solution of the steady flow of a viscous incompressible fluid past a fixed circular cylinder. The solutions obtained using this method produce some interesting results. For instance, an analytic verification of the critical Reynolds number 𝑅 𝑑 for which a standing vortex first appears behind the cylinder is given for the first time and found to be 𝑅 𝑑 2 . 4 . Since these values of the critical Reynolds number are beyond the range of validity of Oseen theory, no analytic verification of this value had previously been given. As a check on the accuracy of the solutions, the calculated drag coefficients at 6th-order approximation are found to agree reasonably well with experimental measurements for 𝑅 𝑑 3 0 which is considerably larger than the 𝑅 𝑑 1 results currently available using other analytic techniques. This buttresses the usefulness of the homotopy analysis method (HAM) as an important tool in solving highly nonlinear problems.

1. Introduction

The viscous incompressible flow past a circular cylinder is a canonical problem which has many practical application areas. It forms a basis for the understanding of bluff body flow which is of significant interest to the research community. The application of a new analytic technique to the cylinder problem is appropriate for the following reasons. First, it forms a benchmark for the application of new mathematical methods. Also various estimates exists in the literature for the critical value of 𝑅 𝑑 for which a standing vortex pair first appears behind the cylinder. Since these values of the Reynolds number lie beyond the range of validity of Oseen theory, obtaining a new method valid for a much higher Reynolds number could aid in providing theoretical confirmations of the actual critical value. Finally, many numerical techniques have been successfully applied in finding approximate solutions to this problem [1] however its full understanding still remains the subject of analytical research. Different analytic studies have been carried out for this problem. These include the far field expansion of the Oseen equations by Filon as well as Shankar [2], the series truncation method of Mandujano and Peralta-Fabi [3], perturbation techniques [4], the Adomians decomposition method [5], the 𝛿 -expansion method [6], and Lyapunov artificial small parameter method [7]. All these approaches produce results which are in agreement with the results of Lamb [8]. Due to the singular nature of the problem, however, no solutions to the general equations have yet been found as regular perturbation techniques becomes singular even for 𝑅 𝑑 1 . A relatively new analytic approach based on the homotopy analysis method is proposed to solve this problem. Different from perturbation techniques, the homotopy analysis method does not depend upon any small or large parameters and thus is valid for most nonlinear problems in science and engineering. Besides, it logically embraces other nonperturbation techniques such as Lyapunov small parameter method, the 𝛿 -expansion method, and Adomians decomposition method. In this paper, we explore the flexibility of the method and choose our boundary conditions in the deformation process appropriately in order that Stokes paradox can be circumvented; this enables us to obtain results which are valid for Reynolds numbers that lie beyond the range of validity of Oseen theory. The governing Navier-Stokes equations expressed in terms of a stream function are first reduced to an infinite set of fourth-order linear partial differential equations using homotopy analysis [923]. The equations are then transformed into a set of linear ordinary differential equations in the radial variable by expanding the flow variables as an infinite series of an elementary transcendental function. This is solved in closed form using the method of variation of parameters. As a check on the accuracy of the solutions obtained, the drag coefficients at 6 th-order approximation are calculated. The calculated drag coefficients are found to agree reasonably well with published experimental as well as numerical results for Reynolds number as high as 30. Furthermore, for higher Reynolds number, the theoretical estimate of the pressure drag tends to become constant while the frictional drag decreases proportionately to the square root of the Reynolds number. Also with regard to the detailed flow patterns our calculations show that for Reynolds number 𝑅 𝑑 below 2.4, a nonseparated flow takes place past the cylinder; this critical value at which the standing vortex pair first appears behind the cylinder is consistent with existing experimental results.

2. Theoretical Framework

Suppose we wish to solve a system of nonlinear partial differential equations given by 𝜓 = 0 ( 2 . 1 ) subject to the boundary conditions 𝜓 ( 𝑎 ) = 𝜓 𝑎 , 𝜓 ( 𝑏 ) = 𝜓 𝑏 ( 2 . 2 ) for some 𝑎 , 𝑏 and 𝜓 𝑛 . Define a homotopy or embedding given by Ψ = 𝑛 × [ ] 1 , 0 𝑛 ( 2 . 3 ) such that Ψ ( 𝑥 , 0 ) = 𝜑 ( 𝑥 ) , Ψ ( 𝑥 , 1 ) = 𝜓 ( 𝑥 ) , ( 2 . 4 ) where 𝜑 ( 𝑥 ) is a solution to the linear problem Δ 0 𝜑 ( 𝑥 ) = 0 ( 2 . 5 ) which can be easily calculated and therefore acts as our starting point or initial guess. Ψ is assumed to be smooth and acts on the solution space by continuously deforming 𝜑 (the starting solution) into 𝜓 the final solution through an embedding or homotopy given by Ψ ( 𝑥 , 𝜆 ) = ( 1 𝜆 ) 𝜑 ( 𝑥 ) + 𝜆 𝜓 ( 𝑥 ) . ( 2 . 6 ) Basically one attempts to trace an implicitly defined curve 𝑐 ( 𝑠 ) in the parametric set of this homotopy map from 𝜑 to a solution point 𝜓 . In other to ensure that this curve intersects the homotopy level, 𝜆 = 1 at finite points, it is sufficient to impose suitable boundary conditions which essentially prevent the curve from running to infinity before intersecting the homotopy level 𝜆 = 1 or from returning to the initial point. One imposes on Ψ the conditions Ψ ( 𝑎 , 𝜆 ) = 𝜑 𝑎 𝜓 + 𝜆 𝑎 𝜑 𝑎 , Ψ ( 𝑏 , 𝜆 ) = 𝜑 𝑏 𝜓 + 𝜆 𝑏 𝜑 𝑏 , ( 2 . 7 ) where 𝜑 ( 𝑎 ) = 𝜑 𝑎 and 𝜑 ( 𝑏 ) = 𝜑 𝑏 . In view of (2.6), (2.7) there exists a family of partial differential operators Δ ( 𝑥 , 𝜆 ) with Δ ( 𝑥 , 0 ) = Δ 0 ( 𝑥 ) , Δ ( 𝑥 , 1 ) = Δ ( 𝑥 ) ( 2 . 8 ) defined by the embedding ( 1 𝜆 ) Δ 0 [ ] ( 𝑥 ) + 𝜆 Δ ( 𝑥 ) , 𝜆 0 , 1 . ( 2 . 9 ) Thus the corresponding homotopy equation is Δ ( 𝑥 , 𝜆 ) Ψ ( 𝑥 , 𝜆 ) = 0 ( 2 . 1 0 ) or by using (2.6) Δ 0 ( 1 𝜆 ) Ψ ( 𝑥 , 𝜆 ) + 𝜆 Δ Ψ ( 𝑥 , 𝜆 ) = 0 . ( 2 . 1 1 ) Equation (2.11) is referred to as the zeroth-order deformation equation. Equation (2.6) defines an implicit relationship between 𝜑 and 𝜓 . In order to make this relationship explicit we assume that the functions Ψ ( 𝑥 , 𝜆 ) are smooth enough about 𝜆 = 0 and expand in a Taylor series about 𝜆 = 0 in the interval [ 0 , 1 ] . Thus we write Ψ ( 𝑥 , 𝜆 ) = Ψ ( 𝑥 , 0 ) + 𝑚 = 1 Ω [ 𝑚 ] 𝜆 𝑚 ! 𝑚 , ( 2 . 1 2 ) where Ω [ 𝑚 ] = 𝜕 𝑚 Ψ ( 𝑥 , 𝜆 ) 𝜕 𝜆 𝑚 | | | | 𝜆 = 0 , 𝑚 0 , ( 2 . 1 3 ) using (2.4) and substituting 𝜓 𝑚 Ω ( 𝑥 ) = [ 𝑚 ] ( 𝑥 ) 𝑚 ! ( 2 . 1 4 ) we have Ψ ( 𝑥 , 𝜆 ) = 𝜙 ( 𝑥 ) + 𝑚 = 1 𝜓 𝑚 ( 𝑥 ) 𝜆 𝑚 . ( 2 . 1 5 ) Assuming the above Taylor series is convergent at 𝜆 = 1 , we have Ψ ( 𝑥 , 1 ) = 𝜙 ( 𝑥 ) + 𝑚 = 1 𝜓 𝑚 ( 𝑥 ) ( 2 . 1 6 ) or 𝜓 ( 𝑥 ) = 𝜙 ( 𝑥 ) + 𝑚 = 1 𝜓 𝑚 ( 𝑥 ) . ( 2 . 1 7 ) The defining equations for 𝜓 𝑚 can be obtained by differentiating the zeroth-order deformation equations 𝑚 times with respect to 𝜆 then setting 𝜆 = 0 and finally dividing by 𝑚 ! . At the boundary points we have 𝜓 ( 𝑎 ) = 𝜙 ( 𝑎 ) + 𝑚 = 1 Ψ 𝑚 ( 𝑎 ) ( 2 . 1 8 ) or 𝑚 = 1 Ψ 𝑚 ( 𝑎 ) = 𝜙 ( 𝑎 ) 𝜓 ( 𝑎 ) ( 2 . 1 9 ) similarly 𝑚 = 1 Ψ 𝑚 ( 𝑏 ) = 𝜙 ( 𝑏 ) 𝜓 ( 𝑏 ) . ( 2 . 2 0 )

3. Problem Definition

For steady two-dimensional flow past a circular cylinder, the governing Navier Stokes equations can be written in terms of the stream function 𝜓 ( 𝑟 , 𝜃 ) and in dimensionless form as 4 𝜓 = R e 𝑟 𝜕 𝜓 𝜕 𝜕 𝜃 𝜕 𝑟 𝜕 𝜓 𝜕 𝜕 𝑟 𝜕 𝜃 2 𝜓 , ( 3 . 1 ) where R e = 𝑎 𝑈 / 𝜈 is the Reynolds number, 𝑎 the radius of the cylinder, 𝑈 the uniform stream velocity at infinity, and 𝜈 the kinematic viscosity. All variables in the above equations are nondimensional. The boundary conditions are given by the no slip condition on the surface of the cylinder 𝜓 ( 1 , 𝜃 ) = 𝜕 𝜓 ( 1 , 𝜃 ) 𝜕 𝑟 = 0 ( 3 . 2 ) and the uniform stream condition at infinity l i m 𝑟 𝜕 𝜓 ( 𝑟 , 𝜃 ) 𝜕 𝑟 = s i n 𝜃 , l i m 𝑟 1 𝑟 𝜕 𝜓 ( 𝑟 , 𝜃 ) 𝜕 𝜃 = c o s 𝜃 . ( 3 . 3 ) In order to obtain an analytic drag formula, the nonlinear partial differential equations (3.1) to (3.3) are mapped into a set of linear subsystem by means of the homotopy technique; these are then solved in closed form using the method of variation of parameters with the aid of the symbolic computation software MATHCAD. Under Stokes approximation, (3.1) reduces to the biharmonic equation 4 𝜓 = 0 . ( 3 . 4 ) Both Stokes and Lamb have pointed out that if the attempt is made to find the steady motion produced by the translation of a cylinder with constant velocity through a liquid of infinite extent on the basis of (3.4), it is impossible to satisfy all the conditions (3.2) and (3.3). It seems that no steady “creeping flow” of a viscous fluid past an infinite cylinder is possible [19]. As a result of this difficulty we choose our boundary conditions in the deformation process appropriately in order that Stokes paradox can be circumvented. A suitable solution to (3.4) is [20] 𝐴 𝜓 = 𝑟 + 𝐵 𝑟 + 𝐶 𝑟 l n ( 𝑟 ) + 𝐷 𝑟 3 s i n 𝜃 , ( 3 . 5 ) where 𝐴 , 𝐵 , 𝐶 , 𝐷 are arbitrary constants to be chosen appropriately in order that the boundary conditions (3.2) and (3.3) are satisfied.

4. Method of Solution

Defining a nonlinear Navier Stokes operator Δ as Δ 𝜓 4 𝜓 R e 𝑟 𝜕 𝜓 𝜕 𝜕 𝜃 𝜕 𝑟 𝜕 𝜓 𝜕 𝜕 𝑟 𝜕 𝜃 2 𝜓 ( 4 . 1 ) and an arbitrary linear operator Δ 0 to be Δ 0 𝜓 4 𝜓 . ( 4 . 2 ) Then by using (4.1) and (4.2) in (2.11), the zeroth-order deformation equations can be written as ( 1 𝜆 ) Δ 0 Ψ ( 𝑟 , 𝜃 , 𝜆 ) 𝜆 Δ Ψ ( 𝑟 , 𝜃 , 𝜆 ) = 0 , ( 4 . 3 ) where 𝑟 1 , 𝜆 [ 0 , 1 ] , 𝜋 𝜃 𝜋 , and < 0 . is an auxiliary parameter often used in HAM to accelerate convergence of the series solutions. In this paper we will simply set = 1 and use the so-called homotopy-Padé technique to get better approximations for our solutions. The associated boundary conditions are given by Ψ ( 1 , 𝜃 , 𝜆 ) = 𝜕 Ψ ( 𝑟 , 𝜃 , 𝜆 ) 𝜕 𝑟 𝑟 = 1 = 0 , l i m 𝑟 Ψ ( 𝑟 , 𝜃 , 𝜆 ) = s i n 𝜃 , l i m 𝑟 Ψ ( 𝑟 , 𝜃 , 𝜆 ) 𝑟 = c o s 𝜃 . ( 4 . 4 ) Differentiating the resulting zeroth order deformation equations 𝑚 times with respect to 𝜆 then setting 𝜆 = 0 and finally dividing by 𝑚 ! , we obtain the 𝑚 th-order deformation equations given as 4 Ψ 𝑚 𝜒 ( 𝑟 , 𝜃 ) = 𝑚 1 4 Ψ 𝑚 1 ( 𝑟 , 𝜃 ) + R e 𝑟 𝑚 1 𝑙 = 0 𝜕 Ψ 𝑙 𝜕 𝜕 𝜃 𝜕 𝑟 𝜕 Ψ 𝑙 𝜕 𝜕 𝑟 𝜕 𝜃 2 Ψ 𝑚 𝑙 1 ( 𝑟 , 𝜃 ) ( 4 . 5 ) with the boundary conditions Ψ 𝑚 ( 1 , 𝜃 ) = 𝜕 Ψ 𝑚 ( 𝑟 , 𝜃 ) 𝜕 𝑟 𝑟 = 1 = 0 , l i m 𝑟 Ψ 𝑚 ( 𝑟 , 𝜃 ) = 0 , l i m 𝑟 Ψ 𝑚 ( 𝑟 , 𝜃 ) 𝑟 = 0 , 𝑚 1 , ( 4 . 6 ) where the coefficient 𝜒 𝑚 is defined by 𝜒 𝑚 = 0 , i f 𝑚 1 , 1 , i f 𝑚 2 . ( 4 . 7 ) We observe that (4.5) gives a system of linear partial differential equations which can be solved analytically for each 𝑚 . The deformations are given by Ψ ( 𝑟 , 𝜃 , 𝜆 ) = Ψ ( 𝑟 , 𝜃 , 0 ) + 𝑚 = 1 𝜓 [ 𝑚 ] 0 𝜆 𝑚 ! 𝑚 , ( 4 . 8 ) where 𝜓 [ 𝑚 ] 0 = 𝜕 𝑚 Ψ ( 𝑥 , 𝜃 , 𝜆 ) 𝜕 𝜆 𝑚 𝜆 = 0 , 𝑚 0 , ( 4 . 9 ) substituting 𝜓 𝑚 𝜓 ( 𝑟 , 𝜃 ) = [ 𝑚 ] 0 ( 𝑟 , 𝜃 ) 𝑚 ! ( 4 . 1 0 ) we have Ψ ( 𝑟 , 𝜃 , 𝜆 ) = 𝜓 0 ( 𝑟 , 𝜃 ) + 𝑚 = 1 𝜓 𝑚 ( 𝑟 , 𝜃 ) 𝜆 𝑚 . ( 4 . 1 1 ) Assuming the above Taylor series is convergent at 𝜆 = 1 , we have 𝜓 ( 𝑟 , 𝜃 ) = 𝜓 0 ( 𝑟 , 𝜃 ) + 𝑚 = 1 𝜓 𝑚 ( 𝑟 , 𝜃 ) . ( 4 . 1 2 ) The above expression gives a direct relationship between the initial approximation 𝜓 0 ( 𝑟 , 𝜃 ) and the final solution 𝜓 ( 𝑟 , 𝜃 ) through the unknown terms 𝜓 𝑚 ( 𝑟 , 𝜃 ) ( 𝑚 = 1 , 2 , 3 , ) , whose governing equations are given by (4.5). In view of (3.5) we set 𝜓 0 𝐴 ( 𝑟 , 𝜃 ) = 𝑟 + 𝐵 𝑟 + 𝐶 𝑟 l n ( 𝑟 ) + 𝐷 𝑟 3 s i n 𝜃 . ( 4 . 1 3 ) In view of the boundary conditions and the form of the initial guess (4.13) as well as the symmetry of the flow, it is reasonable to assume that 𝜓 𝑚 ( 𝑟 , 𝜃 ) = 𝑘 = 1 𝑓 𝑚 , 𝑘 ( 𝑟 ) s i n 𝑘 𝜃 . ( 4 . 1 4 ) Substituting (4.14) in (4.5) we deduce the governing equations 𝑓 𝑖 𝑣 𝑚 , 𝑛 2 ( 𝑟 ) + 𝑟 𝑓 𝑖 𝑖 𝑖 𝑚 , 𝑛 ( 𝑟 ) 2 𝑛 2 + 1 𝑟 2 𝑓 𝑖 𝑖 𝑚 , 𝑛 ( 𝑟 ) + 2 𝑛 2 + 1 𝑟 3 𝑓 𝑖 𝑚 , 𝑛 𝑛 ( 𝑟 ) + 2 𝑛 2 4 𝑟 4 𝑓 𝑚 , 𝑛 ( 𝑟 ) = 𝑞 𝑚 , 𝑛 ( 𝑟 ) 𝑟 1 , 𝑚 1 , ( 4 . 1 5 ) where 𝑞 𝑚 , 𝑛 is given by 𝜒 𝑚 𝑓 1 𝑖 𝑣 𝑚 1 , 𝑛 2 ( 𝑟 ) + 𝑟 𝑓 𝑖 𝑖 𝑖 𝑚 1 , 𝑛 ( 𝑟 ) 2 𝑛 2 + 1 𝑟 2 𝑓 𝑖 𝑖 𝑚 1 , 𝑛 ( 𝑟 ) + 2 𝑛 2 + 1 𝑟 3 × 𝑓 𝑖 𝑚 1 , 𝑛 𝑛 ( 𝑟 ) + 2 𝑛 2 4 𝑟 4 𝑓 𝑚 1 , 𝑛 ( 𝑟 ) R e 4 𝜋 𝑘 = 1 𝑗 = 1 𝑚 1 𝑙 = 0 𝑗 𝑓 𝑖 𝑙 , 𝑘 ( 𝑓 𝑟 ) 𝑖 𝑖 𝑚 1 𝑙 , 𝑗 ( 1 𝑟 ) + 𝑟 𝑓 𝑖 𝑚 1 𝑙 , 𝑗 ( 𝑗 𝑟 ) 2 𝑟 𝑓 𝑚 1 𝑙 , 𝑗 ( × 𝑟 ) 𝜋 𝜋 + s i n 𝑛 𝜃 s i n 𝑘 𝜃 c o s 𝑗 𝜃 𝑑 𝜃 𝑘 = 1 𝑗 = 1 𝑚 1 𝑙 = 0 𝑘 𝑓 𝑙 , 𝑘 𝑓 ( 𝑟 ) 𝑖 𝑖 𝑖 𝑚 1 𝑙 , 𝑗 1 ( 𝑟 ) + 𝑟 𝑓 𝑖 𝑖 𝑚 1 𝑙 , 𝑗 𝑗 ( 𝑟 ) 2 + 1 𝑟 𝑓 𝑖 𝑚 1 𝑙 , 𝑗 + 𝑗 ( 𝑟 ) 2 𝑟 2 𝑓 𝑚 1 𝑙 , 𝑗 ( 𝑟 ) 𝜋 𝜋 . s i n 𝑛 𝜃 s i n 𝑗 𝜃 c o s 𝑘 𝜃 𝑑 𝜃 ( 4 . 1 6 ) Note that 𝑓 𝑖 𝑣 𝑚 , 𝑗 ( 𝑟 ) in (4.15) and (4.16) denotes the 4 th-order derivative of 𝑓 𝑚 , 𝑗 ( 𝑟 ) with respect to 𝑟 . In view of (4.13) and (4.14) we set 𝑓 0 , 1 𝐴 ( 𝑟 ) = 𝑟 𝐷 + 𝐵 𝑟 + 𝐶 𝑟 l n ( 𝑟 ) + 𝑟 3 . ( 4 . 1 7 ) We observe that (4.15) is a well-known ordinary differential equation known as Euler differential equation, which is much easier to solve than the partial differential equation (4.5). The corresponding homogeneous Euler equation 𝑟 4 𝐹 𝑖 𝑣 𝑚 , 𝑛 ( 𝑟 ) + 2 𝑟 3 𝐹 𝑖 𝑖 𝑖 𝑚 , 𝑛 ( 𝑟 ) 2 𝑛 2 𝑟 + 1 2 𝐹 𝑖 𝑖 𝑚 , 𝑛 ( 𝑟 ) + 2 𝑛 2 + 1 𝑟 𝐹 𝑖 𝑚 , 𝑛 ( 𝑟 ) + 𝑛 2 𝑛 2 𝐹 4 𝑚 , 𝑛 ( 𝑟 ) = 0 , 𝑟 1 , 𝑚 1 , ( 4 . 1 8 ) has the general solution 𝐹 𝑚 , 𝑛 ( 𝑟 ) = 𝐶 1 𝑟 𝑛 + 𝐶 2 𝑟 𝑛 + 𝐶 3 𝑟 ( 𝑛 2 ) + 𝐶 4 𝑟 ( 𝑛 + 2 ) , ( 4 . 1 9 ) where 𝐶 1 , 𝐶 2 , 𝐶 3 , 𝐶 4 are constants. The general solution to (4.15) is then given as 𝑓 𝑚 , 𝑛 ( 𝑟 ) = 𝐶 1 𝑟 𝑛 + 𝐶 2 𝑟 𝑛 + 𝐶 3 𝑟 ( 𝑛 2 ) + 𝐶 4 𝑟 ( 𝑛 + 2 ) + 𝑓 𝑚 , 𝑛 ( 𝑟 ) , ( 4 . 2 0 ) where 𝑓 𝑚 , 𝑛 ( 𝑟 ) is the complimentary function. In order to obtain expressions for 𝑓 𝑚 , 𝑛 ( 𝑟 ) , we use the method of variation of parameters. We assume for each 𝑚 and each 𝑛 1 that the complimentary function 𝑓 𝑚 , 𝑛 ( 𝑟 ) is of the form 𝑟 𝑛 𝑓 ( 𝑟 ) + 𝑟 𝑛 𝑔 ( 𝑟 ) + 𝑟 ( 𝑛 2 ) ( 𝑟 ) + 𝑟 ( 𝑛 + 2 ) 𝑤 ( 𝑟 ) , ( 4 . 2 1 ) where 1 𝑓 ( 𝑟 ) = 8 𝑞 e x p ( l n ( 𝑟 ) 𝑛 ) 𝑚 , 𝑛 [ ] 𝑞 ( 1 + 𝑛 ) 𝑟 𝑛 𝑑 𝑟 , 𝑔 ( 𝑟 ) = 𝑚 , 𝑛 1 8 e x p ( l n ( 𝑟 ) 𝑛 ) 𝑛 ( 𝑛 1 ) 𝑟 𝑑 𝑟 , ( 𝑟 ) = 8 𝑞 𝑚 , 𝑛 e x p ( l n ( 𝑟 ) 𝑛 ) 𝑛 ( 𝑛 1 ) 𝑟 3 𝑤 𝑞 𝑑 𝑟 , ( 𝑟 ) = 𝑚 , 𝑛 8 e x p ( l n ( 𝑟 ) 𝑛 ) ( 1 + 𝑛 ) 𝑛 𝑟 3 𝑑 𝑟 . ( 4 . 2 2 ) Note that when 𝑛 = 1 the solutions for arbitrary 𝑛 given by (4.19) are not linearly independent. A linearly independent solution to (4.18) for 𝑛 = 1 is given as 𝐹 1 , 1 𝐴 ( 𝑟 ) = 𝑟 𝐷 + 𝐵 𝑟 + 𝐶 𝑟 l n ( 𝑟 ) + 𝑟 3 . ( 4 . 2 3 ) Thus for 𝑛 = 1 we assume a complimentary function of the form 𝑟 1 𝑓 ( 𝑟 ) + 𝑟 ( 𝑝 ( 𝑟 ) + l n ( 𝑟 ) 𝑤 ( 𝑟 ) ) + 𝑟 3 ( 𝑟 ) , ( 4 . 2 4 ) where 𝑞 𝑓 ( 𝑟 ) = 𝑚 , 1 𝑞 1 6 𝑑 𝑟 , 𝑝 ( 𝑟 ) = 𝑚 , 1 4 𝑟 2 𝑞 𝑑 𝑟 , ( 𝑟 ) = 𝑚 , 1 1 6 𝑟 4 𝑞 𝑑 𝑟 , 𝑤 ( 𝑟 ) = l n ( 𝑟 ) 𝑚 , 1 4 𝑟 2 𝑑 𝑟 . ( 4 . 2 5 ) The initial guess, which is a stokes approximant, does not satisfy the boundary conditions at infinity but contains an unknown coefficient 𝐶 . The value of 𝐶 is determined by enforcing the higher-order approximations together with the initial guess to satisfy the boundary condition at infinity. Hence, we solve the set of differential equations (4.15) one after the other, subject to the following boundary conditions: 𝑓 𝑚 , 𝑛 ( 1 ) = 𝑓 𝑖 𝑚 , 𝑛 ( 1 ) = 0 , 𝑛 1 , 𝑚 0 , ( 4 . 2 6 ) l i m 𝑟 𝑝 𝑚 = 0 𝑓 𝑚 , 1 ( 𝑟 ) 𝑟 = 1 , l i m 𝑟 𝑝 𝑚 = 0 𝑓 𝑖 𝑚 , 1 ( 𝑟 ) = 1 , ( 4 . 2 7 ) l i m 𝑟 𝑓 𝑚 , 𝑛 ( 𝑟 ) 𝑟 = 0 , l i m 𝑟 𝑓 𝑖 𝑚 , 𝑛 ( 𝑟 ) = 0 , 𝑛 1 , 𝑚 0 , ( 4 . 2 8 ) so that 𝐶 is chosen in order to make (4.27) satisfied. The series (4.14) hence (4.27) must be truncated at some fixed value 𝑝 . As 𝑝 is increased it is expected that the solutions converge to the solutions of the original equations (3.1). This has not been proven yet. Clearly the complexity of the differential equations increases rapidly with 𝑝 as such the number of terms we can compute becomes restricted by the ability of the computational software to handle the long terms. We obtain for the first few values of 𝑓 𝑛 , 𝑚 ( 𝑟 ) , where the constants 𝐴 , 𝐵 , 𝐶 , 𝐷 have been appropriately chosen to satisfy the boundary conditions (4.26) to (4.28): 𝑓 0 , 1 𝐶 ( 𝑟 ) = 2 1 𝑟 , 𝑓 𝑟 + 2 𝑟 l n ( 𝑟 ) 1 , 1 𝑓 ( 𝑟 ) = 0 , 1 , 2 1 ( 𝑟 ) = 9 6 𝑟 2 l n ( 𝑟 ) 2 + 1 7 5 7 6 𝑟 2 1 1 1 0 𝑟 7 0 𝑟 3 + 1 1 9 2 𝑟 4 1 l n ( 𝑟 ) 1 3 0 𝑟 7 1 7 3 5 0 𝑟 3 + 1 1 9 2 0 𝑟 6 1 7 7 5 3 5 6 4 4 8 0 𝑟 2 + 1 9 4 6 0 8 𝑟 4 + 2 2 9 1 3 9 𝐶 5 6 4 4 8 0 0 2 𝑓 R e , 2 , 1 ( 𝑟 ) = 9 9 1 9 9 0 4 0 1 1 3 8 1 3 6 5 7 1 2 6 3 0 6 6 9 7 9 9 5 8 7 8 4 0 0 0 0 0 𝑟 1 2 1 2 2 4 6 7 1 0 3 8 3 9 1 8 7 + 1 1 4 5 7 3 8 4 9 7 6 8 7 5 5 2 0 0 0 0 0 𝑟 1 5 3 6 𝑟 5 + 1 1 5 3 6 𝑟 3 1 1 9 2 0 𝑟 4 + 1 1 9 2 0 𝑟 4 l n ( 𝑟 ) 4 + 1 8 0 𝑟 5 4 5 7 5 2 6 4 𝑟 6 3 3 4 9 1 1 8 8 1 6 0 0 𝑟 4 + 6 0 7 5 7 6 0 0 𝑟 2 1 5 7 2 4 5 7 6 𝑟 7 + 4 2 7 3 6 8 6 4 𝑟 5 + 7 9 6 0 𝑟 3 l n ( 𝑟 ) 3 + 2 2 9 1 3 9 2 2 5 7 9 2 0 0 2 7 4 7 5 1 9 4 8 3 2 6 4 0 𝑟 6 + 1 4 9 2 1 5 0 4 𝑟 8 8 2 0 7 1 2 3 + 1 8 6 0 3 3 6 0 0 𝑟 2 3 9 2 3 1 1 0 3 2 1 9 2 0 0 𝑟 7 1 5 7 3 8 4 0 0 𝑟 9 1 0 0 5 6 3 1 2 1 6 7 6 0 3 2 0 𝑟 5 + 5 8 7 1 4 1 3 7 6 3 2 0 0 0 𝑟 2 + 2 8 6 7 8 2 2 9 5 4 1 9 0 0 8 0 0 𝑟 3 1 4 2 8 6 1 0 9 4 7 4 1 6 3 2 0 0 𝑟 4 l n ( 𝑟 ) 2 + 2 5 8 1 2 1 7 5 7 2 1 6 7 6 0 3 2 0 0 0 𝑟 7 4 4 4 7 7 0 9 6 3 2 0 𝑟 1 0 3 1 7 4 2 0 9 8 3 6 7 7 3 7 6 0 0 0 0 𝑟 2 + 6 3 3 2 3 2 3 7 2 7 1 9 9 1 4 8 5 4 4 0 0 0 𝑟 4 2 4 5 1 5 9 2 5 8 0 4 8 0 0 𝑟 9 7 6 4 5 9 9 + 3 6 1 2 6 7 2 0 𝑟 6 1 3 6 8 6 4 0 𝑟 1 1 2 2 8 3 8 0 4 6 1 2 3 8 9 7 8 2 5 2 8 0 𝑟 6 4 5 0 5 6 2 0 8 1 1 3 0 0 5 6 1 9 2 0 0 𝑟 5 3 2 5 4 7 5 4 1 1 6 5 0 2 8 0 9 6 0 0 𝑟 3 + 2 2 9 1 3 9 + 1 1 2 8 9 6 0 0 2 4 5 1 4 8 7 3 2 6 8 4 0 8 9 6 0 0 𝑟 8 l n ( 𝑟 ) + 4 6 7 7 2 7 4 0 9 9 9 1 8 5 8 9 0 9 4 0 9 2 8 0 0 𝑟 8 + 3 5 3 1 7 2 0 3 2 0 𝑟 1 3 + 2 2 9 1 3 9 3 3 8 6 8 8 0 0 1 7 1 9 4 7 6 1 9 3 1 5 2 0 0 𝑟 1 1 + 2 1 5 7 2 9 9 1 5 2 4 0 9 6 0 0 𝑟 3 2 7 8 8 0 9 2 9 1 2 5 4 0 0 0 0 0 𝑟 2 + 2 3 5 2 4 0 3 7 4 9 1 7 7 4 9 7 6 0 0 𝑟 1 0 + 3 3 7 0 7 8 8 5 0 6 3 5 2 0 2 2 4 7 6 8 0 0 0 0 𝑟 7 5 4 1 3 2 8 6 0 4 0 7 5 3 7 7 0 1 0 6 8 8 0 0 0 𝑟 6 + 3 4 4 2 5 4 8 5 4 1 6 1 7 9 8 1 4 4 𝑟 5 1 9 0 7 8 3 8 6 5 6 0 𝑟 1 2 1 2 8 5 5 0 4 7 6 1 9 3 1 5 2 0 0 0 𝑟 9 + 2 9 6 1 8 5 5 2 4 2 1 7 2 0 9 1 0 5 9 7 1 2 0 0 0 0 𝑟 4 𝐶 3 R e 2 + 1 + 1 2 𝑟 2 𝑓 𝑟 𝑟 l n ( 𝑟 ) 𝐶 , 2 , 2 ( 𝑓 𝑟 ) = 0 , 2 , 3 ( 𝑟 ) = 4 2 0 5 3 2 7 4 2 0 3 2 5 4 7 8 1 0 3 2 5 8 0 2 6 5 4 0 1 9 1 5 8 0 1 6 0 0 0 0 𝑟 3 + 1 1 8 3 6 4 0 1 2 7 4 9 6 0 9 9 6 7 + 1 1 1 4 6 7 4 6 5 3 6 2 7 6 2 6 2 9 1 2 0 𝑟 3 1 3 6 2 8 8 𝑟 4 + 3 8 1 9 2 𝑟 5 1 7 5 7 6 0 𝑟 3 l n ( 𝑟 ) 3 + 3 1 6 7 1 4 9 2 1 6 7 6 0 3 2 0 0 𝑟 3 9 8 9 1 2 4 7 4 0 0 𝑟 6 7 2 7 9 5 2 5 6 𝑟 4 + 3 4 9 7 2 2 9 3 7 6 0 𝑟 5 + 8 3 2 7 6 4 8 0 𝑟 7 l n ( 𝑟 ) 2 + 2 0 4 7 1 1 2 1 1 2 0 0 𝑟 8 + 1 5 3 1 1 2 9 9 6 3 3 7 9 2 0 0 𝑟 5 6 8 6 4 8 4 8 3 3 5 5 3 2 4 6 8 4 8 0 0 0 𝑟 6 + 3 1 3 8 7 0 7 2 𝑟 9 + 3 5 9 1 2 9 9 2 8 9 7 2 8 0 0 𝑟 7 1 7 9 5 7 4 3 1 1 9 2 0 3 6 0 9 6 0 𝑟 4 + 4 8 1 6 8 6 9 1 2 6 0 1 1 2 3 8 4 0 0 𝑟 3 + 2 2 9 1 3 9 1 3 1 7 1 2 0 0 0 𝑟 3 + l n ( 𝑟 ) 1 5 1 9 3 6 8 3 1 1 1 5 6 0 5 5 0 4 0 0 0 𝑟 5 + 1 1 7 0 3 7 1 4 4 1 3 7 8 2 4 4 9 1 1 3 6 0 0 0 0 𝑟 6 1 8 3 9 9 8 6 1 7 2 7 6 5 9 5 2 0 0 0 𝑟 2 2 4 3 1 5 5 4 1 3 5 9 1 1 7 1 5 8 4 0 0 𝑟 8 2 8 3 0 0 2 8 3 1 9 5 0 8 4 2 8 8 0 0 0 𝑟 7 + 1 6 7 0 0 3 5 2 0 2 2 4 7 6 8 0 𝑟 9 2 1 2 3 3 1 9 8 1 3 2 4 8 0 0 𝑟 1 0 2 9 0 0 4 4 9 4 4 5 9 2 4 1 9 6 5 4 8 0 9 6 0 0 𝑟 4 1 9 1 2 7 5 2 5 1 2 0 𝑟 1 1 𝐶 3 R e 2 , ( 4 . 2 9 ) and so on. All these values are valid for the entire flow field. Unfortunately we are unable to give approximations higher than order six due to the nonlinearity 𝑞 𝑚 𝑛 ( 𝑟 ) in (4.15) which makes the terms of 𝑓 𝑚 , 𝑘 ( 𝑟 ) increase almost exponentially and the computational software MATHCAD 2000 cannot display the results.

5. Analytic Drag Formula

The force on the cylinder's surface is calculated using 𝐹 = 𝑆 𝜏 𝑟 𝑟 c o s 𝜃 𝜏 𝑟 𝜃 s i n 𝜃 𝑑 𝑠 , ( 5 . 1 ) where the component of the stress tensor is evaluated at the cylinder's surface 𝑆 . Introducing the stream function 𝜓 and rewriting the components in dimensionless form we have that the nondimensional drag coefficient 𝐶 𝐷 can be written as 𝐶 𝐷 2 = R e 𝜋 0 𝜕 3 𝜓 𝜕 𝑟 3 𝜕 + 2 2 𝜓 𝜕 𝑟 2 𝑟 = 1 s i n 𝜃 𝑑 𝜃 , ( 5 . 2 ) which upon substitution of (3.2) in (5.2) can be written as 𝐶 𝐷 𝜋 = R e 𝑚 = 0 𝑓 𝑖 𝑖 𝑖 𝑚 , 1 ( 1 ) + 2 𝑓 𝑖 𝑖 𝑚 , 1 ( 1 ) . ( 5 . 3 ) Making the substitutions we have our 4 th-order drag formula to be 𝐶 𝐷 = 9 2 3 7 0 8 6 4 9 3 9 6 5 5 5 4 5 5 9 6 7 6 2 7 4 3 6 8 0 3 5 0 1 8 9 2 9 4 2 0 6 9 3 5 8 4 6 1 6 0 8 2 2 2 2 7 5 3 9 1 8 9 7 3 4 2 0 5 6 3 4 2 8 8 9 8 0 5 3 0 9 1 7 1 4 4 9 4 2 6 7 8 7 6 1 4 7 3 5 1 0 4 2 3 3 1 4 4 9 5 5 0 6 7 9 5 7 2 4 8 2 5 3 4 2 8 2 3 8 0 0 0 0 0 0 0 0 0 R e 3 𝐶 5 + 6 3 6 9 1 8 9 7 8 1 1 7 1 0 8 6 2 2 4 5 8 1 2 8 5 2 1 9 6 5 5 8 6 9 1 8 7 0 4 9 3 1 2 7 4 7 5 2 0 0 0 0 0 0 R e 𝐶 3 𝜋 , ( 5 . 4 ) where 𝐶 satisfies 6 9 9 9 7 8 5 8 4 7 4 2 9 2 3 7 1 1 0 8 1 2 8 6 2 9 8 5 6 3 2 8 0 0 5 7 0 8 8 6 1 7 9 0 0 4 9 4 9 4 0 7 7 2 9 8 4 0 4 7 9 1 5 8 0 5 0 3 0 6 4 7 6 0 5 5 3 5 8 7