High-Reynolds-number flow past a shear-free circular cylinder

A semi-analytical study of the dynamics of incompressible flow of a Newtonian fluid past a circular cylinder with zero surface shear stress has been performed. Essentially, the modified flow field is considered to be a perturbed flow field over the irrotational flow description around a circular cylinder in an uniform cross flow. At finite Reynolds number, the perturbation over irrotational flow field is confined within a thin boundary layer, rear stagnation and wake regions. Equations governing the flow-field in boundary layer are linearized Navier-Stokes equations which is solved by adopting analytical route. Solution is shown to be superposition of infinite self-similar terms. Behind the circular cylinder, a stagnation region is formed where the perturbation and irrotational flow velocity components are shown to be of equal order. Vorticity produced in boundary layer is convected downstream through stagnation region to form a narrow wake. Similar to boundary layer analysis, governing equation in wake is also linearized to facilitate the analytical study. Solution in wake region is obtained by matching the solution from boundary layer near the stagnation region. Far wake solution is shown to be consistent with the leading order drag coefficient estimated from the momentum deficit in the far wake. Next, drag coefficient is estimated by calculating the net-dissipation over the entire flow field. Asymptotic expansion of drag coefficient in terms of boundary layer thickness is shown to have a logarithmic singularity. We show that at very high Reynolds number, flow past a circular cylinder with shear-free boundary will cause a net increase in dissipation of energy over the viscous irrotational flow, which stands in stark contrast to the case of spherical geometry (Moore[1]). Results from present analysis is validated against computational results.


I. INTRODUCTION
Characterization of the flow past a bluff body, such as a circular cylinder, has always been a problem of great interest to fluid dynamicists. Flow past a circular cylinder has been studied quite extensively in history (Strouhal [2], Von Karman [3], Williamson [4]). This problem often serves as a benchmark to understand principal phenomena in bluff body flows such as flow separation and vortex shedding. Typically, the analysis of flow past a bluff body proceeds with the imposition of a no-slip boundary condition on the surface. Although the application of no-slip condition at solid -liquid interface is now widely accepted, it should be understood that the no-slip is just an approximation that holds under the most of the flow conditions and is not something that can be derived from first principle (Batchelor [5], Leal [6], Lauga et al. [7]). Indeed, there are numerous scenarios, both natural and technological, where no-slip boundary condition at the surface fails to describe the true fluid velocity at the surface of a solid body. Leaves of many plants, one typical example is of lotus leaf, exhibit water repulsion property and provide finite slip to water droplets moving on it's surface (Barthlott and Neinhuis [8], Bhushan and Jung [9]). This property has initially served as a source of inspiration for developing superhydrophobic surfaces. Superhydrophobic surfaces are capable of bringing appreciable drag reduction in both laminar and turbulent flow regimes (Ou et al. [10], Daniello et al. [11], Muralidhar et al. [12]) and can provide intermediate to extremely large slip lengths (Rothstein [13]). Given the tremendous growth in research on superhydrophobic surfaces, there is a need for characterization of flow past a bluff body with prescribed Navier slip boundary condition on its surface, and quantification of reduction in hydrodynamic loads acting on it.
Another, instance where we encounter vanishingly small shear-stress or perfect-slip condition, is the case of multiphase flows of air and water (Harlow et al. [14], Scardovelli and Zaleski [15]). Simplest example of above setting is the motion of a spherical gas bubble through a liquid. Levich [16] was the first to theoretically analyze the drag on a spherical bubble at high Reynolds number. He assumed that the flow is only slightly perturbed from irrotational flow and found the drag on the bubble up to first order through calculation of net dissipation in the flow-field. Moore [1] extended his drag calculations up to one more order by finding the correct velocity field near the surface of bubble through boundary layer analysis and generalized it to oblate ellipsoids (Moore [17]). Bubbles are prone to deformation; therefore Moore's analysis of spherical bubble is not applicable at very high Reynolds number. Vakarelski et al. [18] suggested an alternative to this problem, by creating a robust vapor layer covering the whole surface of bluff body moving in a liquid, using Leidenfrost phenomenon. They got a very thin stable wake behind the sphere even at Re = 2 × 10 5 , as predicted by Moore in the analysis of spherical bubble.
For a typical bluff body flow, the introduction of slip at the interface lead to lower momentum losses; a consequence of enhanced momentum transport in the boundary layer enabled by hydrodynamic slip. This reduced momentum loss is signified by a drop in normal velocity gradients and hence decrease in vorticity production at the interface. This reduction in the vorticity production accompanied by its augmented convection downstream results in decrease of vorticity accumulation behind the bluff body and thereby diminishing the strength of the unsteady vortical structures in wake and recirculation near the rear stagnation point (You and Moin [19], Legendre et al. [20], Seo and Song [21], Li et al. [22]). The enhanced momentum transport also facilitates increased pressure recovery in the downstream side which leads to substantial reduction in pressure drag and as a result significant drop in overall hydrodynamic loads are achieved.
Legendre et al. [20] through numerical solutions showed that for a circular cylinder maximum non-dimensional vorticity at the surface should exceed a critical value of (ω max s a/U ∞ 5) to trigger vortex shedding. Except in the Stokes regime, no theory is generally available for determining the drag experienced by a bluff body with no slip or finite slip boundary. However, in flow over a shear-free cylinder , the maximum non-dimensional vorticity on the surface at high Reynolds number is asymptotic to a value of 4 and therefore the flow will remain laminar and attached at all Reynolds number.
In this paper we perform a semi-analytical study of flow past a shear-free circular cylinder at high Reynolds number. The paper is organized as follows: section II provides the description of the setup of the flow configuration and the governing equations. The theoretical analysis is outlined for the boundary layer, rear stagnation and wake region along with our major findings and validation against numerical simulations are given in section III. Section IV discusses in detail the derivation of drag coefficient from previous analysis and validation against numerical simulation before concluding the work in section V.

II. PROBLEM DEFINITION
We consider the viscous incompressible flow of a Newtonian fluid past a uniform circular cylinder of infinite span and radial dimension a. The interface between the cylinder wall and the fluid is under perfect hydrodynamic slip condition. The problem is best described in the plain-polar coordinates (r, θ). The equations governing the flow are the steady-state, 2-D Navier-Stokes equations: along with the no through-flow and shear-free conditions at the cylinder surface where u = (u r , u θ ) is the velocity field, p denotes the hydrodynamic pressure and, ρ and µ are density and dynamic viscosity of the fluid respectively. The non-dimensional parameter characterizing the flow; Reynolds number is defined based on the cylinder diameter and is given by Re = ρU ∞ D/µ, where U ∞ denotes the free-stream velocity and D = 2a is the diameter of the cylinder. The flow approaches uniform stream at large distances from the cylinder which gives the far-field conditions as follows where p ∞ denotes the free-stream pressure.

III. ANALYSIS
The velocity field and the pressure for the potential flow past a cylinder is given by will satisfy the Navier-Stokes equations (1) and no through-flow (2a) at cylinder surface as well as the conditions at infinity (3). But, this flow-field fails to satisfy the zero shear-stress condition (2b) at the cylinder surface, which at high Reynolds numbers will lead to the development of a thin boundary layer across which the relative change in tangential stress is O(1) as shown by Moore [23] for spherical bubble. The vorticity produced in the boundary layer is convected downstream through a stagnation region resulting in a narrow wake behind the cylinder. In the forthcoming sections we present the theoretical analysis in each of these regions viz. the boundary layer, the stagnation zone and the wake in detail.

III.1. Boundary layer analysis
To analyze the boundary layer let us define a radial coordinate system scaled by the cylinder radius and attached to the cylinder surface: y = (r − a)/a, along with a new azimuthal coordinate with origin at the front stagnation point: φ = π − θ.
Let us write the flow variables as the superposition of potential flow and correction over it, in the newly defined coordinate system where u y and u φ are the radial and azimuthal potential flow velocity components in (y, φ) coordinate system while p denotes the potential flow pressure. The quantitiesũ y ,ũ φ andp denote the corresponding corrections. Use (5b) and (4b) in (2b), and rewrite the shear-free condition at the interface If δ denotes the boundary layer thickness scaled with radius of the cylinder a then within the boundary layer (∂/∂y) varies as O(1/δ). The potential flow velocity components u φ and u y in the boundary layer are O(1) and O(δ) respectively. Using the above facts it can be easily deduced from (6) that in the boundary layerũ φ is O(δ). Further applying continuity equation we getũ y to be O(δ 2 ). For deriving the boundary layer equations let us define the coordinate system (y * , φ) with appropriate stretching in the radial direction, here y * = y/δ. Also, we will non-dimensionalise the total velocity components as follows: The azimuthal momentum equation in terms of (7) and coordinates (y * , φ) is: In the boundary layer there exists a balance between convective terms and diffusive terms, and hence the quantity µ/ρU ∞ aδ 2 = 2/(Reδ 2 ) should be O(1). Without loss of generality let us choose δ = 2/Re.
The radial momentum equation in terms of (7) and coordinates (y * , φ) is: From equation (9) one can see that in the boundary layer the potential flow pressure (p) provides the required centripetal force for the flow to turn around the cylinder (see the underlined terms) and therefore leaving pressure perturbationsp to be O(δ 2 ).
Using (5) in (8) along with the potential flow variables (4a-c) and retaining only the leading order terms (boundary layer approximations), we obtain the governing equation of non-dimensionalised correction velocity in azimuthal The boundary conditions forũ * φ are obtained from (2b) and (3) with only leading order terms retained as Note that the condition (11a) follows from the symmetry of the flow.
The solution to (10) subject to boundary conditions (11a-c) may be expressed as an infinite superposition of selfsimilar terms i.e.ũ * φ = ∞ n=0 w nũ * φ,n , where w n 's are constants. Eachũ * φ,n is a self-similar solution of the equation (10) which satisfy the condition ∂ũ * φ,n /∂y * = f n (φ) at y * = 0 instead of (11c). It turns out that one possible form for f n (φ) is f n (φ) = (sin(φ/2)) 2n+1 where n = 0 or any positive integer. Similarity solutionũ * φ,n with this f n (φ) which also satisfies the first two boundary conditions (11a,b) is where η = √ 2y * cos φ 2 . The right hand side of original boundary condition for the radial derivative ofũ * φ (11c) can be expressed as a linear combination of f n (φ) and the linearity of the system implies that the same superposition of u * φ,n will solve (10) subject to (11a-c) which leads to the following final solution forũ * The above solution is in the form of superposition of infinite self-similar solutions as opposed to a single self-similar solution in case of spherical geometry (Moore [1]).
Vorticity of O(1) gets generated in the boundary layer (ω bl ) and has the following non-dimensional form.

III.2. Rear stagnation region analysis
As the flow approaches the rear stagnation point, it takes a sharp turn and the gradients in azimuthal direction become comparable to the normal gradients which leads to breakdown of boundary layer approximation. We adopt the (y, θ) coordinates (y as defined in §III.1) in this region so that the origin of our local reference frame coincides with the rear stagnation point. Boundary layer thickens as it approaches the stagnation region as a result of decrease in vorticity gradient. This leads to convection dominated dynamics over viscous diffusion near the stagnation region. Harper [24] analyzed the viscous layers in which vorticity changes appreciably but not the velocity (similar to our analysis in boundary layer). This analysis showed that the region where viscous terms remain negligible in comparison of convective terms extends up to a distance from the stagnation region into the boundary layer which is of O (1) . He also showed that the boundary layer assumptions are valid only beyond a distance from rear stagnation point which is where the boundary layer assumptions are valid and the region is inviscid. Vorticity in the final stage of boundary layer can be calculated as follows.
In the final stage of the boundary layer, where the distortion in stream-function (ψ) is still negligible in comparison of potential flow stream-function (ψ), the net stream-function (ψ) which is the sum of above two, is asymptotic to ψ = ψ +ψ ≈ ψ ≈ 2U ∞ ayθ. Using this asymptotic form, η can be re-written as where ψ * s = ψ/U ∞ aδ is the non-dimensional form of stream-function near stagnation region. This expression of η, (16) along with (15) gives the expression of vorticity in the final stage of boundary layer. Because this region is inviscid, meaning vorticity is constant along streamline and hence, the same expression will give the vorticity in stagnation region: ω * s is the non-dimensional vorticity in the stagnation region given as ω * s = ω s a/U ∞ . Analyzing the stagnation region is nothing but analyzing the flow in a corner with angle π/2, means both the directions, y and θ are equally important or in other words y ∼ θ. Also, the flow turns around the corner which implies i.e. ψ ∼ yθ. Combining this with fact that ψ is O(Re −1/2 ) (see 16) in stagnation region, gives y, θ ∼ δ 1/2 or in other words size of stagnation region is O(Re −1/4 ). Define the stretched coordinate system (y * s , θ * s ) where y * s = y/δ 1 2 and θ * s = θ/δ 1 2 . Absence of vortex stretching mechanisms suggests that the vorticity level in stagnation zone should be same as that in the boundary layer i.e. ω s ∼ O(1), which can be seen from (17) as well. The velocity corrections should undergo a change of O(Re −1/4 ) across this region to carry the vorticity of O(1); which implies the correction velocity components are O(Re −1/4 ). From (4a,b) it can be shown that in this region the potential velocity components are also of O(Re −1/4 ), which means that correction and potential flow components are of comparable magnitude. Writing the vorticity as a curl of velocity vector and retaining only the significant terms based on scalings in stagnation region, we get u * θ (s) =ũ θ /U ∞ δ 1/2 andũ * y (s) =ũ r /U ∞ δ 1/2 are the non-dimensional correction velocity components in azimuthal and radial directions respectively. We write correction velocity components in terms of distortion in stream-function as followsũ * Use (17) and (19) in (18), we obtain the governing equation of distortion in stream-function as follows ψ * s =ψ/U ∞ aδ and ψ * s = ψ/U ∞ aδ = 2y * s θ * s are the non-dimensional distortion and potential flow stream-function respectively, in the stagnation region. We obtain an elliptic equation with non-linear source term because in the case of two-dimensional planar geometry perturbations are comparable to potential flow in the stagnation region (Harper [24]); whereas this is not the situation in the case of two-dimensional axi-symmetric geometry (Moore [1]), because the magnitude of vorticity decreases due to contraction of vortex-lines as flow approaches stagnation region and therefore perturbations remains smaller in comparison potential flow (Harper and Moore [25]).
The correction flow become unidirectional in the final stage of boundary layer and in the beginning of the wake. This fact along with the no through-flow at cylinder surface and the symmetry condition along the rear axis of symmetry (θ = 0) translates into the following boundary conditions Non-linear equation (20) along with boundary conditions (21a-c) is solved numerically to obtain the flow-field in the stagnation region (details are given in appendix). Fig.1 depicts the tangential velocity correction (ũ φ ) scaled with Re/2 on the surface of cylinder, obtained using analytical expression (13) and the result from stagnation region analysis (19) compared against numerical results. The numerical results were obtained via solving the Navier-Stokes equation in plain polar coordinates using 10 th order compact finite difference schemes (Shukla et al. [26], Shukla and Zhong [27]) along radial coordinate and Fourier spectral schemes in azimuthal coordinate. Temporal discretization is done by a 2 nd order semi-implicit projection scheme (Hugues and Randriamampianina [28]) which uses Adam's Bashforth and Backward differentiation formula (AB/BDI) for nonlinear convective and linear viscous terms respectively. It can be clearly seen from Fig.1 that the numerical results approaches results from boundary layer and stagnation region analysis as Reynolds number increases.

III.3. Wake analysis
The vorticity generated in the boundary layer is convected downstream in a narrow wake behind the cylinder. A right-handed Cartesian coordinate system (x, z) in which either coordinates are scaled with cylinder radius a is chosen for analyzing the flow in this region. The origin of the coordinate system is at the rear stagnation point with x axis pointing in stream-wise direction.
Vorticity level is O(1) same as boundary layer and stagnation region, in the initial stage of wake because of the absence of vortex stretching mechanism. Now if δ w denotes the wake thickness scaled with cylinder radius a, then the correction velocity in stream-wise directionũ x should be O(δ w ) to carry the vorticity produced in boundary layer. Conservation of mass principle gives the magnitude of other component of correction velocity i.eũ z as O(δ 2 w ). The magnitude of potential flow velocity components u x and u z in the wake are easily found to be O(1) and O(δ w ) respectively. Now again by balancing the convective terms with diffusion terms in wake similar to the boundary layer, it can be shown δ w ∼ δ, so we have chosen δ w = δ.
Based on above information let us define appropriate scaled spatial coordinates, potential flow velocity components and correction velocity components as follows: Linearized governing equation for the vorticity in wake region with only the leading order terms retained is given by where u * x0 and u * z0 are the zeroth terms in Taylor series expansion of u * x and u * z about z * = 0 respectively, while ω * w = ω w a/U ∞ is the non-dimensionalised vorticity in the wake. Symmetry of the flow about z * = 0 along with irrotationality condition at far field leads to the following boundary conditions In addition to that solution of ω * w should be such that it would give rise to a constant momentum deficit in far wake. The general solution of (23) which satisfies (24a,b) along with the far wake momentum deficit constraint is Here λ n and c n are constants whose values are to be obtained by comparing the asymptotic forms of ω * w with ω * bl in the rear stagnation region. Now the vorticity in the beginning of wake is equal to the vorticity in final stage of boundary layer because the flow in stagnation region is symmetric in y and θ. Matching of general solution of vorticity in wake (25) with the vorticity solution (15) gives the final expression for ω * w as follows: As we go far downstream,i.e. as x → ∞, (26) reduces to Away from the cylinder surface the wake solution exhibits the correct decay rate that is consistent with the leading order drag coefficient which is 16π/Re. Further the above expression can be shown to be in agreement with the vorticity in the far wake of a general 2-D body as given by (Schlichting and Gersten [29]).

IV. VISCOUS DISSIPATION AND DRAG COEFFICIENT
An expression for the net drag experienced by the circular cylinder is obtained by computing the total viscous dissipation over the entire flow field. This method is based on mechanical energy balance. It does not require pressure and stress distribution over the surface of the cylinder. Calculation of drag by force balance on the cylinder requires correction in pressure in the boundary layer which is of O(Re −1 ) and hard to find. Also, this method gives the information only about leading order term in drag coefficient (Kang and Leal [30]).
The non-dimensional total viscous dissipation (effective drag coefficient C D ) for flow past a circular cylinder with a prescribed tangential surface velocity u θ (a, θ) is given by the following expression (Shukla and Arakeri [31]) Calculation of the first integral in (28) is straight-forward. We have evaluated this integral upto O(Re −3/2 ) which is sufficient to calculate the leading order correction in C D over the same at infinite Reynolds number limit.
The first integral on the right hand side in (29) corresponds to the leading order term in drag coefficient. Using the expression of (u θ ) and (ũ φ ) from (4b and 13) and calculating second integral numerically on the right hand side of (29), we get Due to the symmetry of flow, the second integral in (28) may be written as where Ω u represents the region (r, θ) : 0 ≤ θ ≤ π, a ≤ r < ∞. Next, we will show that vorticity square integral (31) near stagnation region exhibits a logarithmic singularity. Therefore, using the singularity subtraction technique we write 2µ ρU 3 ∞ a Ωu In (32), Ω s denotes stagnation region in the upper half of the flow. The first integral (I) is calculated over the whole domain using only the analytical solutions from boundary layer (14) and wake (26). In this integral, we are actually taking into account only the asymptotic values of the vorticity from the boundary layer or the wake in the stagnation region and not the actual vorticity in the stagnation region itself. Therefore, the second integral (II) which is free of logarithmic singularity and calculated only near the stagnation region, accounts for the difference.
For the calculation of first integral (I), we have two solutions; one from the boundary layer (14) and other is from the wake (26) region. But there is no certain boundary to divide domain into two regions of validity for these two solutions. However, we know that these two solutions have same asymptotic form in the rear stagnation region. Therefore, it is plausible to assume that the curve dividing the domains passes through the rear stagnation point. Finally, it can be shown that for any curve C which originates from the rear stagnation point to divide the domain Ω u into two regions Ω bl and Ω w , (I) will exhibit the same leading order term. Using this idea, let us split the double integral (31) as shown below: To calculate the integral as in (33), we only need the information about the initial profile of the curve C since outside the thin stagnation region, vorticity quickly drops to zero. Assume that the asymptotic equation of curve C near the stagnation region is θ ∼ Ay b . Curve C after originating from the stagnation point should proceed in upright direction because we are considering only the domain above the axis of symmetry and the wake solution is valid only in downstream portion of rear stagnation point. It poses constraints that A > 0 and b > 1/2. Also, when b > 1/2, the equation for C in (x, z) coordinates is z ∼ Ax b . Using these asymptotic formulas of curve C, we compute the  integral (33) up to the leading order terms as Note that above result is independent of the constants A and b and hence shape of the dividing curve C. Now the integral II which is calculated only near the stagnation region can be written as follows. 2µ Since, we do not have an analytical solution in the stagnation region, above integral can only be obtained numerically to get We could have constructed a single uniformly valid solution of vorticity using solution from the boundary layer, stagnation and wake to calculate (31). But, as explained by Van Dyke [32], if the solutions in different regions share common regions of validity, then the leading order terms will remain independent of method of constructing the composite solution. Finally, adding (30), (34) and (36) to get the expression for drag coefficient as: Define correction term over dissipation from potential flow with appropriate scaling as  (37) and numerical calculation is less than the drag from viscous potential flow analysis, but both results are asymptotic to it as Reynolds number increases. Fig. 2b shows the comparison of correction in drag coefficient as defined in (38) from the analytical expression (39) against computational results. Same trend is obtained in both the results, that the net increase in dissipation over the potential flow goes from negative to positive at certain Reynolds number and also, numerical results are asymptotic to analytical results at high Reynolds number. Both the results show, increase in net dissipation over potential flow at very large Reynolds number. This result is significantly different from the case of shear-free sphere (Moore [1]), where the net dissipation is always lower than the potential flow even at high Reynolds number.

V. CONCLUSIONS
In this study, we have developed a theoretical framework for the flow over a shear-free circular cylinder. At large Reynolds numbers, the flow under consideration is only a slightly perturbed version of viscous irrotational flow. A very thin boundary layer develops over the cylinder surface across which the corrections over irrotational flow decays to zero. The corrections within the boundary layer is found to be a linear combination of self-similar terms. Behind the cylinder, a stagnation region is formed where the corrections and base flow components are of comparable order of magnitude, which is not the situation in case of axis-symmetric bluff bodies. The vorticity generated in boundary layer gets convected in a narrow wake through the rear stagnation region. A full analysis of wake region has been carried out and the solution obtained was found to possess a constant momentum deficit in far wake, which corresponds to the leading order drag coefficient, 16π/Re. The drag coefficient (C D ) was obtained by estimating the dissipation of energy over the entire flow field. In stark contrast to the case of spherical geometry, for which the drag reduction produced by shear-free flow was greater compared to the viscous irrotational flow at all Reynolds number values, in the present case the same holds true only for low and moderate Reynolds numbers. At large Reynolds numbers, contrary to our expectation, drag is lower for the viscous irrotational flow than the shear-free flow.
Appendix A: Solution of equation (10) First we show that the following equation along with boundary conditionsũ * φ,n = 0 at φ = 0, exhibits a similarity solution ∀ n 0 and n ∈ Z. Assume thatũ * φ,n = H n (φ)F n (η), where η = y * /g(φ). Inserting this form into (A1), we obtain Similarity requires the terms in the brackets to be some constants. Without loss of generality we have chosen The only solution of (A4) such that the final solution of (A1) will satisfy the boundary condition (A2a) is Boundary condition (A2c) requires c n = 2n + 2 and d n = Finally the governing equation of F n is obtained by inserting (A4) into (A3) and the boundary condition (A2b) will transform into Solution of (A7) which satisfies the boundary condition (A8) is given by Above can be shown using mathematical induction. First we show that for n = −1, (A9) satisfies (A7).
Which indeed is the solution of (A7) for n = −1. (A10) is a simple derivation of Craig's formula (Craig [33]) which another form of Gaussian distribution. Next we show that following relationship holds between F n+1 (η) and F n (η) Integrating the equation (A7) twice one can obtain Using relationship (A11) into (A12) shows that if (A9) is the solution of (A7) for some integer n then it also holds for n + 1 which then completes the proof. Also, We cannot directly use the limit η = 0 in the integrand because the integrand is singular at α = π. Therefore, we first simplify above using integration by parts and then we use the limit η = 0 in the resultant well-behaved integrand here B(x, y) represent beta function. Our original boundary condition (11c) can be written as where Because of the linearity of the equation (10), we superpose the similarity solutionsũ * φ,n using the weights (A15b) to obtain the solution of scaled azimuthal correction velocity component (ũ * φ ) as follows where η is given by Appendix B: Numerical solution of equation (20) Consider the following transformation Write equation (20) in (σ, τ ) coordinates is We have to solve above equation with following boundary conditions Boundary condition (B3a) means that cylinder wall is impermeable, (B3b) means that outside the thin stagnation layer distortion in stream function become constant, (B3c) means that flow is symmetric about the line θ * s = y * s and (B3d) means that the flow in the end stage of boundary layer is unidirectional.
To solve (B2) we use iterative process. We start with a guess valueψ * s = 0 everywhere, then in each iteration we use the known value ofψ * s on R.H.S. of (B2) from previous iteration and solve forψ * s by treating as an unknown on L.H.S. of (B2) along with boundary conditions (B3a-d). We keep doing it untilψ * s converges. To solve the Poisson's equation that we obtain in each iteration, we use 2nd order finite difference method on exponentially stretched grids in both σ and τ directions. The resultant system of linear equations is then solved using Gaussian Elimination method by taking the advantage of sparsity of obtained matrix.
Tangential velocity at the surface of cylinder is given bỹ Derivative in above equation at σ = 0 is obtained using 6th order forward finite difference method.
Using converged value ofψ * s , vorticity in stagnation region is obtained as follows Correction in vorticity square integral near stagnation region (35) in (σ, τ ) coordinates can be written as Above integral is calculated using Simpson's 2-D method.
Appendix C: Calculation of integrals in equation (34) First find out the leading order terms in in the limit of δ → 0 + . Given conditions are also, M 0 , N 0 , M ′ 0 and N ′ 0 are bounded. Performing the integration w.r.t. t in (C1) will give Divide the integral into two parts. Integration limits in part one is from s = 0 to s = ǫ and in second part it is from s = ǫ to s = s f . Choose such value of ǫ that δ 1 n+1 ≪ ǫ and ǫ ≪ 1. Because of such order of ǫ, in part one integrand can be simplified by using the limit ǫ → 0 and at the same time in second part integrand can be simplified by taking the limit ǫ n+1 /δ → ∞. Simplification of part one is as follows From (A10) we have the following relation Use this relation in (C4) and perform the integration w.r.t. ξ, we obtain E 1 is the exponential integral function. From Abramowitz and Stegun [34] we have, E 1 (x) ∼ −γ − ln x as x → 0, where γ = 0.577 is Euler's constant. Using this, simplify the limits in (C6) Using the fact that ǫ n+1 ≫ δ and, erf(x) ∼ 1 − exp(−x 2 )/ √ πx and E 1 (x) ∼ exp(−x)/x as x → ∞, we get the leading order terms in (C7) as follows: ln ǫ in the leading order term of I 1 is the reason of order of error term in second step of (C4). Now simplify I 2 Integrand in above integration can be expanded using the fact that s n+1 /δ ≫ 1 Use the singularity subtraction technique to simplify (C10) as follows where s 0 s f is some O(1) number can be chosen based on convenience. Now, the integrand in second integral is bounded as s → 0, therefore the leading order terms in I 2 is given as Finally adding (C7) and (C12) we get the leading order terms in (C1) One can see that in above expression δ has explicitly separated out, and the integrals which are independent of δ can be calculated numerically if required.
Next will show that the integrals in (34) can be written using the form (C1) and therefore can be calculated in the similar manner as described above. Square of vorticity from boundary layer in θ coordinate can be written as Using this one can write the first integral in (34) as (C16) Also, the second integral in (34) can be written as Now to calculate vorticity square integral over whole domain we will again use the same technique as for the case of cylinder. Choose curve C as φ ∼ A((r−a)/a) b . Similar to case cylinder, constraints are A > 0 and b > 1/2. If b > 1/2, equation of curve C in (x, m) coordinates is m/a ∼ A((x− a)/a) b or in (x * , s * ) coordinates is s * ∼ 3/δA(x * /δ) b+1/2 . Using this curve C calculate the vorticity square integral. 2µ y=0 ω 2 bl sin φ dy dφ In this case, as the vorticity in stagnation region is given by asymptotic value of vorticity from the boundary layer.
Integral in (D12) is calculated numerically to obtain II = 89.47 One thing to notice here is that the leading order terms in both the integral (I) and (II) are independent of shape of curve C, this is because contribution of vorticity square integral near the stagnation region is negligible (O(Re −11/6 )) (Moore [1]). Also, if we use the axisymmetric far wake equation instead of equation (D3) but solve with same boundary conditions (D4), solution for vorticity then obtained will give the same leading order term in II as obtained in (D13). This is why Moore was able to use the far wake equation directly to simplify integral (4.27) in his paper to obtain drag coefficient. Moore's approach cannot be used in case of cylinder because the vorticity solution in wake, obtained by considering transport due to potential flow rather than uniform flow differ in the leading order term. Finally, adding (D9), (D11) and (D13) will reproduce the Moore's result of drag on a shear-free sphere