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

Authors: E. O. Ifidon
Journal: Journal of Applied Mathematics
 (2009)

---

## 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 numberRd for which a standing vortex first appears behind the cylinder is given for the first time and found to be Rd≼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 Rd≃30 which is considerably larger than the Rd≃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.


---

## Body

## 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 ofRd 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 Rd≪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 [9–23]. 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 6th-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 Rd 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(2.1)∇ψ=0
subject to the boundary conditions
(2.2)ψ(a)=ψa,ψ(b)=ψb
for some a,b∈ℝ and ψ∈ℝn. Define a homotopy or embedding given by
(2.3)Ψ=ℝn×[1,0]→ℝn
such that
(2.4)Ψ(x,0)=φ(x),Ψ(x,1)=ψ(x),
where φ(x) is a solution to the linear problem
(2.5)Δ0φ(x)=0
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
(2.6)Ψ(x,λ)=(1-λ)φ(x)+λψ(x).
Basically one attempts to trace an implicitly defined curve c(s) 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)Ψ(a,λ)=φa+λ(ψa-φa),Ψ(b,λ)=φb+λ(ψb-φb),
where φ(a)=φa and φ(b)=φb. In view of (2.6), (2.7) there exists a family of partial differential operators Δ(x,λ) with
(2.8)Δ(x,0)=Δ0(x),Δ(x,1)=Δ(x)
defined by the embedding
(2.9)(1-λ)Δ0(x)+λΔ(x),λ∈[0,1].
Thus the corresponding homotopy equation is
(2.10)Δ(x,λ)Ψ(x,λ)=0
or by using (2.6)
(2.11)Δ0(1-λ)Ψ(x,λ)+λΔΨ(x,λ)=0.
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 Ψ(x,λ) are smooth enough about λ=0 and expand in a Taylor series about λ=0 in the interval [0,1]. Thus we write
(2.12)Ψ(x,λ)=Ψ(x,0)+∑m=1∞Ω[m]m!λm,
where
(2.13)Ω[m]=∂mΨ(x,λ)∂λm|λ=0,m≥0,
using (2.4) and substituting
(2.14)ψm(x)=Ω[m](x)m!
we have
(2.15)Ψ(x,λ)=ϕ(x)+∑m=1∞ψm(x)λm.
Assuming the above Taylor series is convergent at λ=1, we have
(2.16)Ψ(x,1)=ϕ(x)+∑m=1∞ψm(x)
or
(2.17)ψ(x)=ϕ(x)+∑m=1∞ψm(x).
The defining equations for ψm can be obtained by differentiating the zeroth-order deformation equations m times with respect to λ then setting λ=0 and finally dividing by m!. At the boundary points we have
(2.18)ψ(a)=ϕ(a)+∑m=1∞Ψm(a)
or
(2.19)∑m=1∞Ψm(a)=ϕ(a)-ψ(a)
similarly
(2.20)∑m=1∞Ψm(b)=ϕ(b)-ψ(b).


## 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ψ(r,θ) and in dimensionless form as
(3.1)∇4ψ=Rer(∂ψ∂θ∂∂r-∂ψ∂r∂∂θ)∇2ψ,
where Re=aU∞/ν is the Reynolds number, a the radius of the cylinder, U∞ 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
(3.2)ψ(1,θ)=∂ψ(1,θ)∂r=0
and the uniform stream condition at infinity
(3.3)limr→∞∂ψ(r,θ)∂r=sinθ,limr→∞1r∂ψ(r,θ)∂θ=cosθ.
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
(3.4)∇4ψ=0.
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]
(3.5)ψ=[Ar+Br+Crln(r)+Dr3]sinθ,
where A,B,C,D 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.1)Δψ≡∇4ψ-Rer(∂ψ∂θ∂∂r-∂ψ∂r∂∂θ)∇2ψ
and an arbitrary linear operator Δ0 to be
(4.2)Δ0ψ≡∇4ψ.
Then by using (4.1) and (4.2) in (2.11), the zeroth-order deformation equations can be written as
(4.3)(1-λ)Δ0Ψ(r,θ,λ)-λℏΔΨ(r,θ,λ)=0,
where r≥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
(4.4)Ψ(1,θ,λ)=∂Ψ(r,θ,λ)∂r∥r=1=0,limr→∞Ψ(r,θ,λ)=sinθ,limr→∞Ψ(r,θ,λ)r=cosθ.
Differentiating the resulting zeroth order deformation equations m times with respect to λ then setting λ=0 and finally dividing by m!, we obtain the mth-order deformation equations given as
(4.5)∇4Ψm(r,θ)=(χm-1)∇4Ψm-1(r,θ)+Rer∑l=0m-1(∂Ψl∂θ∂∂r-∂Ψl∂r∂∂θ)∇2Ψm-l-1(r,θ)
with the boundary conditions
(4.6)Ψm(1,θ)=∂Ψm(r,θ)∂r∥r=1=0,limr→∞Ψm(r,θ)=0,limr→∞Ψm(r,θ)r=0,m≥1,
where the coefficient χm is defined by
(4.7)χm={0,ifm≼1,1,ifm≽2.
We observe that (4.5) gives a system of linear partial differential equations which can be solved analytically for each m. The deformations are given by
(4.8)Ψ(r,θ,λ)=Ψ(r,θ,0)+∑m=1∞ψ0[m]m!λm,
where
(4.9)ψ0[m]=∂mΨ(x,θ,λ)∂λm∥λ=0,m≥0,
substituting
(4.10)ψm(r,θ)=ψ0[m](r,θ)m!
we have
(4.11)Ψ(r,θ,λ)=ψ0(r,θ)+∑m=1∞ψm(r,θ)λm.
Assuming the above Taylor series is convergent at λ=1, we have
(4.12)ψ(r,θ)=ψ0(r,θ)+∑m=1∞ψm(r,θ).
The above expression gives a direct relationship between the initial approximation ψ0(r,θ) and the final solution ψ(r,θ) through the unknown terms ψm(r,θ)(m=1,2,3,…), whose governing equations are given by (4.5). In view of (3.5) we set
(4.13)ψ0(r,θ)=[Ar+Br+Crln(r)+Dr3]sinθ.
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
(4.14)ψm(r,θ)=∑k=1∞fm,k(r)sinkθ.
Substituting (4.14) in (4.5) we deduce the governing equations
(4.15)fm,niv(r)+2rfm,niii(r)-(2n2+1)r2fm,nii(r)+(2n2+1)r3fm,ni(r)+n2(n2-4)r4fm,n(r)=qm,n(r)r≽1,m≽1,
where qm,n is given by
(4.16)(χm-1)(fm-1,niv(r)+2rfm-1,niii(r)-(2n2+1)r2fm-1,nii(r)+(2n2+1)r3×fm-1,ni(r)+n2(n2-4)r4fm-1,n(r))-Re4π{∑k=1∞∑j=1∞∑l=0m-1(jfl,ki(r)[fm-1-l,jii(r)+1rfm-1-l,ji(r)-j2rfm-1-l,j(r)]×∫-ππsinnθsinkθcosjθdθ)+∑k=1∞∑j=1∞∑l=0m-1(kfl,k(r)[fm-1-l,jiii(r)+1rfm-1-l,jii(r)-(j2+1)rfm-1-l,ji(r)]+j2r2fm-1-l,j(r)∫-ππsinnθsinjθcoskθdθ)}.
Note that fm,jiv(r) in (4.15) and (4.16) denotes the 4th-order derivative of fm,j(r) with respect to r. In view of (4.13) and (4.14) we set
(4.17)f0,1(r)=[Ar+Br+Crln(r)+Dr3].
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.18)r4Fm,niv(r)+2r3Fm,niii(r)-(2n2+1)r2Fm,nii(r)+(2n2+1)rFm,ni(r)+n2(n2-4)Fm,n(r)=0,r≽1,m≽1,
has the general solution
(4.19)Fm,n(r)=C1r-n+C2rn+C3r-(n-2)+C4r(n+2),
where C1, C2, C3, C4 are constants. The general solution to (4.15) is then given as
(4.20)fm,n(r)=C1r-n+C2rn+C3r-(n-2)+C4r(n+2)+fm,n*(r),
where fm,n*(r) is the complimentary function. In order to obtain expressions for fm,n*(r), we use the method of variation of parameters. We assume for each m and each n≻1 that the complimentary function fm,n*(r) is of the form
(4.21)r-nf(r)+rng(r)+r-(n-2)h(r)+r(n+2)w(r),
where
(4.22)f(r)=∫(-18exp(ln(r)n)qm,n[(1+n)rn])dr,g(r)=∫(-qm,n8exp(ln(r)n)n(n-1)r)dr,h(r)=∫(18qm,nexp(ln(r)n)[n(n-1)r3])dr,w(r)=∫(qm,n8exp(ln(r)n)(1+n)nr3)dr.
Note that when n=1 the solutions for arbitrary n given by (4.19) are not linearly independent. A linearly independent solution to (4.18) for n=1 is given as
(4.23)F1,1(r)=[Ar+Br+Crln(r)+Dr3].
Thus for n=1 we assume a complimentary function of the form
(4.24)r-1f(r)+r(p(r)+ln(r)w(r))+r3h(r),
where
(4.25)f(r)=∫(-qm,116)dr,p(r)=∫(-qm,14r2)dr,h(r)=∫(qm,116r4)dr,w(r)=∫(ln(r)qm,14r2)dr.
The initial guess, which is a stokes approximant, does not satisfy the boundary conditions at infinity but contains an unknown coefficient C. The value of C 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:
(4.26)fm,n(1)=fm,ni(1)=0,n≽1,m≽0,(4.27)limr→∞(∑m=0pfm,1(r)r)=1,limr→∞(∑m=0pfm,1i(r))=1,(4.28)limr→∞fm,n(r)r=0,limr→∞fm,ni(r)=0,n≻1,m≻0,
so that C is chosen in order to make (4.27) satisfied. The series (4.14) hence (4.27) must be truncated at some fixed value p. As p 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 p 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 fn,m(r), where the constants A,B,C,D have been appropriately chosen to satisfy the boundary conditions (4.26) to (4.28):
(4.29)f0,1(r)=C2[1r-r+2rln(r)],f1,1(r)=0,f1,2(r)=[196r2ln(r)2+(17576r2-110r-170r3+1192r4)ln(r)-1130r-717350r3+11920r6-17753564480r2+194608r4+2291395644800]C2Re,f2,1(r)=[-99199040113813657126306697995878400000r-1212246710383918714573849768755200000r+(-11536r5+11536r3-11920r4+11920r4)ln(r)4+(-180r-54575264r6-334911881600r4+60757600r2-15724576r7+42736864r5+7960r3)ln(r)3+(22913922579200-27475194832640r6+14921504r8-8207123186033600r+239231103219200r7-15738400r9-1005631216760320r5+58714137632000r2+28678229541900800r3-14286109474163200r4)ln(r)2+(25812175721676032000r7-44477096320r10-3174209836773760000r2+6332323727199148544000r4-24515925804800r9-76459936126720r+61368640r11-22838046123897825280r6-45056208113005619200r53254754116502809600r3+22913911289600+245148732684089600r8)ln(r)+4677274099918589094092800r8+3531720320r13+22913933868800-171947619315200r11+2157299152409600r3-27880929125400000r2+2352403749177497600r10+337078850635202247680000r7-541328604075377010688000r6+344254854161798144r5-19078386560r12-128550476193152000r9+29618552421720910597120000r4]C3Re2+[-12r+12r-rln(r)]C,f2,2(r)=0,f2,3(r)=[4205327420325478103258026540191580160000r3+118364012749609967111467465362762629120r+(-3136288r4+38192r5-175760r3)ln(r)3+(31671492167603200r3-9891247400r6-72795256r4+34972293760r5+83276480r7)ln(r)2+(-204711211200r8+1531129963379200r5-686484833553246848000r6+31387072r9+359129928972800r7-179574311920360960r4+4816869126011238400r3+229139131712000r3]ln(r)+151936831115605504000r5+1170371441378244911360000r6-1839986172765952000r2-24315541359117158400r8-28300283195084288000r7+1670035202247680r9-212331981324800r10-290044944592419654809600r4-19127525120r11]C3Re2,
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 qmn(r) in (4.15) which makes the terms of fm,k(r) 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(5.1)F=∫S(τrrcosθ-τrθsinθ)ds,
where the component of the stress tensor is evaluated at the cylinder's surface S. Introducing the stream function ψ and rewriting the components in dimensionless form we have that the nondimensional drag coefficient CD can be written as
(5.2)CD=-2Re∫0π[(∂3ψ∂r3)+2(∂2ψ∂r2)]r=1sinθdθ,
which upon substitution of (3.2) in (5.2) can be written as
(5.3)CD=-πRe∑m=0∞[fm,1iii(1)+2fm,1ii(1)].
Making the substitutions we have our 4th-order drag formula to be
(5.4)CD=[923708649396555455967627436803501892942069358461608222275391897342056342889805309171449426787614735104233144955067957248253428238000000000Re3C5+636918978117108622458128521965586918704931274752000000ReC3]π,
where C∈ℝ satisfies
(5.5)(69997858474292371108128629856328005708861790049494077298404791580503064760553587328225320193849468139486207709426951494803419822347736747792741280000000000)Re3C5+6260245091238567649156865691931173837409862549504000000Re2C3=1.
The 6th-order drag formula is
(5.6)CD=[(423745417981864617679078623498332027352271339112799555522385048910179519037142905060556090687762380465063908600305420854795501674929953713777012378346150588892986396910104753042147912472232298833897190248711096168716941830424053346724259858599880008971177532410062642147184181594683733597181152000000000000000000)Re5C7+(4756166159482024826680129597541831414741204992000000C6-129868071106241849159270986310108749606119088676708035115110261417811423927048731255659992288578223241045368658233493999845306260143019080130062168658233493999845376000000000000C5)Re3-(85057115567484061097733608448204721773228208884734785068106215827308312803116126133021861629330166984749992845279769617614900325462133114496843408998400000000000)C2Re2-(850571155674840610977473478506810621582730831803116126133021861656014080341202955380348141665474213294936269150054243688852416140293301669847568166400000000000C2+394402210300943309098195591965586918704931274752000000C3)Re+4756166159482024826680129896312747122111807488000000Re]π,
where C∈ℝ satisfies
(5.7)(751742209773879079006241780329833140325659174593679011541205827612190564359971625333892731724070838830929210361897965981390871816360708365195937066880000000000)Re4C6+6238328565020662246476748791931173837409862549504000000Re2C3-(4388745484608464157563997262795351410780492219925295889631044228136107711422285743036333654412657428279038345160183252512877301004951412544355158457293981430085293014030337507220156841921953136738889826814922665770123016509825443200803455591515999280053827065194460375888657601074479784029912000000000000000000)Re6C7=1.
Liao and Cheung proposed the so-called homotopy-Padé method which is more efficient than the traditional Padé approximant [12]. Better results for the drag coefficients may be obtained using this technique. From (4.11) we have the (2m)th approximation
(5.8)∂2Ψ(r,θ,λ)∂r2∥r=1≈∑p=02m∂2ψp(r,θ)∂r2∥r=1λp.
Employing the traditional Padé technique to the power series of the embedding parameter λ, we have the [m,m] Padé approximant
(5.9)∂2Ψ(r,θ,λ)∂r2∥r=1≈1+∑p=12mγm,pλp1+∑p=12mγm,m+pλp,
where γm,p is a coefficient. Setting λ=1 in the above series and using (4.11) we obtain
(5.10)∂2ψ(r,θ)∂r2∥r=1≈1+∑p=1mγm,p1+∑p=1m‍γm,m+p.
Note that in the homotopy-Padé method, one first analyzes the power series in the embedding parameter λ and then sets λ=1. This is the distinction between the traditional Padé approximant and the homotopy-Padé ones.


## 6. Results and Discussion

Figure1 shows a comparison of our 4th- and 6th-order HAM drag formula with previous theoretical [24, 25] as well as experimental [26] results. We observe that our 6th order drag formula agrees well with experimental results in the region Rd≃30 (where Rd=2Re) which is much higher than the Rd≪1 range of previous analytic results [3]. We also observe from Figure 1 that the higher the order of approximation, the better the agreement between our drag formula and experimental data. Thus our drag formula seems to be valid in a larger region of Reynolds number as the number of approximations increases. As noted in [15], this kind of tendency of convergence is very important, especially when a rigorous mathematical proof of convergence cannot be provided. Also we observe from Figure 1 that the [3,3] homotopy-Padé approximate increases further the accuracy of our solutions and can give a remarkably good representation of the function up to Rd≈100 if we could get high enough approximations greater than 6th order. With regards to the detailed flow pattern, we consider the question of the formation of the standing eddie. The equation of the separated streamlines is given by Ψ=0, and the angular coordinate of the point of separation is found by setting r=1 in (4.14) for Ψ. Phenomenologically, the flow presents different regimes depending on the value of the Reynolds number. As the Reynolds number is increased from some initial small value, the lack of fore-and-aft symmetry slowly becomes apparent until a recirculating region develops on the wake, and attached to the cylinder. As Re is increased, a stage is reached at which the separation starts to occur, the angle of separation at which it first occurs will be given by cosθ=1, in which case the fluid starts to separate from the rear generator of the cylinder. The critical value of Re will then be that value which makes
(6.1)∑n=1∞∑m=0∞‍nfm,nii(1)=0.
Owing to the limitations of computer storage space it was found necessary to restrict the number of terms in the series to n=5 and m=6. Our calculations show that for Reynolds number Rd 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. Also Figure 2 suggests that as the Reynolds number increases the pressure drag becomes constant while the frictional drag tends to zero. This behaviour of the drag coefficients is also consistent with experimental observations.Figure 1
Comparison with previous theoretical and experimental results.Figure 2
Drag coefficient as a function ofRd.



---

*Source: 524307-2009-07-01.xml*