---
abstract: 'This paper proposes a solution to Stokes’ paradox for asymptotically uniform viscous flow around a cylinder. The existence of a [*global*]{} stream function satisfying a perturbative form of the two-dimensional NavierStokes equations for low Reynolds number is established. This stream function satisfies the appropriate boundary conditions on both the cylinder and at infinity, but nevertheless agrees with Stokes’ original results at finite radius as the Reynolds number tends to zero. The NavierStokes equations are satisfied to a power-log power of the Reynolds number. The drag on the cylinder is calculated from first principles and the free parameter of the approach can be chosen to give good agreement with data on drag. In this revised working paper we put our approach on a firmer mathematical basis using the Helmholtz-Laplace equation as a linear approximation to the NavierStokes system. In so doing we demonstrate the instability of the original paradox. We also demonstrate the absence of a paradox of Stokes-Whitehead class, and give further theoretical constraints on the free parameters of the model.'
author:
- 'William T. Shaw[^1]'
title: 'A simple resolution of Stokes’ paradox?[^2]'
---

Key Words: Stokes Paradox, Fluid dynamics, Stokes flow, Stream function, Biharmonic equation, Helmholtz equation, Low Reynolds number

Introduction
============

The difficulty in establishing a sensible global solution to the problem of low $R_e$ (Reynolds number) viscous flow around simple objects, where the flow is uniform at infinity, has fascinated applied mathematicians for just over 150 years. Stokes (1851) established that there was no solution to the two-dimensional, steady, incompressible, equations for asymptotically uniform flow around a cylinder in the biharmonic limit. This situation is now routinely described in modern tutorial discussions. See, for example, Chapter 7 of Acheson (1990) for an exercise on Stokes’ paradox and a discussion of the corresponding situation for the sphere. The biharmonic equation can be solved in a neighbourhood of the cylinder (i.e. a circle) by any stream function of the form $$\psi = C\sin\theta(r \log r - \frac{1}{2}r + \frac{1}{2r})$$ but there is no choice of $C$ for which $\psi \sim r \sin\theta$ for large $r$. The analysis of this problem has lead to several classic papers (Oseen, 1910; Lamb, 1911) and its understanding through the use of matched asymptotic expansions (MAE) is one of the triumphs of perturbation theory. The reader is referred to Van Dyke (1964) for his classic survey of the work of Kaplun (1957), Proudman and Pearson (1957) and other key references.

However, the MAE approach, despite its immense power and diversity of expanding applications, does not give a clean resolution of the original difficulty in that such methods rely essentially on computing and then matching solutions to the problem defined on two regions: close to the cylinder and far from the cylinder. The purpose of this paper is to address the problem of finding [*global*]{} solutions for the low Reynolds number limit, i.e. to resolve the original paradox. There are, of course, other approaches to the paradox. Recently, Villas Boas [@mario] has considered the problem from a three-dimensional perspective and points out that there is then no paradox.

Viscous incompressible flow in 2D
=================================

A large class of fluids can be characterized by their density, $\rho$, a scalar field not presumed to be constant, and their dynamic viscosity $\mu$. The flow is characterized by a velocity vector field $\underline{v}$, and an associated scalar pressure field $p$. Conservation of mass is expressed by the continuity equation $$\frac{\partial \rho}{\partial t} + \del.(\rho \vvec) = 0$$ and the conservation of momentum is expressed by the Navier-Stokes equations[^3] $$\rho(\frac{\partial \vvec}{\partial t} + \vvec.\del \vvec) = -\del p + \mu \delsq \vvec$$ If the fluid is incompressible in the sense that $\rho$ is a constant in both time and space, we have the condition: $$\del.\vvec = 0$$ To analyze matters further, we introduce the vorticity vector $$\omvec = \del \times \vvec$$ In the following discussion we demand incompressibility but allow for non-zero vorticity. Using simple identities from vector calculus the Navier-Stokes equations may then be recast in the form $$\rho(\frac{\partial \vvec}{\partial t} - \vvec \times \omvec) + \del (p + \frac{1}{2} \rho \vvec^2
)= -\mu \del \times \omega .$$ Taking the curl of this, we arrive at the vorticity equation $$\frac{\partial \omvec}{\partial t} + \vvec.\del \omvec - \omvec.\del \vvec = \nu \delsq \omvec$$ where the kinematic viscosity $\nu = \mu/\rho$.

The stream function
-------------------

Since the velocity field is divergence-free, we may introduce a vector potential $\Psivec$ such that $$\vvec = \del \times \Psivec$$ and furthermore we may choose it so that it is divergence free: $$\del.\Psivec = 0$$ The vector potential can be reduced to a single function when there is an appropriate symmetry. The resulting object is a stream function. For example, planar 2D flow is obtained by setting (and note that this automatically satisfies satisfies the divergence condition) $$\Psivec = - \Psi(x, y,t) \underline{e}_z$$ Next we note that under the assumption that $\Psi$ satisfies $\del.\Psivec = 0$ $$\omvec = -\delsq \Psivec =  \underline{e}_z \delsq \Psi$$ and the vorticity equation becomes

$$\nu \nabla^4 \Psi - \frac{\partial\ }{\partial t}\nabla^2 \Psi=  \frac{\partial(\Psi, \delsq \Psi)}{\partial(X, Y)}$$

For problems where it is possible to identify a natural length scale $L$ and a natural speed $U$, it is normal practice to perform a non-dimensionalization of the variables and introduce the Reynolds number $R_e = UL/\nu$. Then the NavierStokes equation becomes (after rescaling the variables suitably to $\psi, x, y, \tau$): $$\nabla^4 \psi - R_e  \frac{\partial\ }{\partial \tau}\nabla^2 \psi= R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)}$$ Throughout this paper we shall work in units in which the [*radius*]{} $a$ of the cylinder is taken to be unity. In such units the Reynolds number is here based on the radius[^4] and is given by $R_e = U/\nu$. The two stream functions are related by $\Psi = U \psi$ and furthermore $X=x, Y=y$ etc.

The old historical approach to the limiting case when $R_e \rightarrow 0$ and the flow is time-independent is to take the view that the non-linearities may then be ignored (provided the non-linear term is well behaved) and the time-independent NavierStokes equations reduce to $$\nabla^4 \psi = 0$$ which is the biharmonic limit, also known as Stokes flow. It is now well known (see for example, Chapter 8 of Van Dyke (1964)) that the neglect of the non-linear terms can lead to inconsistencies, as is evidenced by the lack of any solution for asymptotically uniform two-dimensional flow past a cylinder.

Here we introduce the scalar vorticity function $\omega = \nabla^2 \psi$ and write the NavierStokes as the pair $$\begin{split}
\omega & = \nabla^2 \psi\\
\nabla^2 \omega - R_e  \frac{\partial \omega}{\partial \tau}&= R_e \frac{\partial(\psi, \omega)}{\partial(x, y)}
\end{split}$$ If we had exponential growth in vorticity, $\omega = e^{\mu \tau}\hat{\omega}(x,y)$, then the linearized form would be $$\begin{split}
 e^{\mu \tau}\hat{\omega}(x,y) & = \nabla^2 \psi\\
\nabla^2\hat{\omega} - R_e \mu \hat{\omega}&= 0\end{split}$$ or in terms of a single condition: $$\nabla^4 \psi - R_e \mu \nabla^2 \psi=0\ ,$$ which is the Laplace-Helmholtz equation.

The instability of Stokes’ paradox
----------------------------------

Let us focus temporarily on the linearized time-dependent case based on the Laplace-Helmholtz equation. We let $\varepsilon^2 = R_e \mu$. It is an elementary exercise (and we shall give equivalent details in a different context later) to establish that for all $\varepsilon >0$, the Laplace-Helmholtz equation $$\nabla^4 \psi - \varepsilon^2 \nabla^2 \psi=0$$ has a solution of the form $$\psi = e^{\mu \tau}\left[r
-\biggl(1+
\frac{2
   K_1(\varepsilon )}{\varepsilon  K_0(\varepsilon
   )}
\biggr) 
\frac{1}{r}
+\frac{2 K_1(r \varepsilon )}{\varepsilon 
 K_0(\varepsilon )}
 \right] \sin (\theta ) \label{tdsf}$$ that satisfies $$\begin{split}
\psi = 0 = \frac{\partial \psi}{\partial r} \ {\rm if }\  r= 1\\
\psi \sim  e^{\mu \tau} r \sin(\theta) \  {\rm as}\  r \rightarrow \infty
\end{split}$$ In other words, no matter how small the value of $\mu$, the paradox does not exist.

To summarize: If we have some time-dependence in the simplified form of exponential growth in vorticity the linearization is precisely the Laplace-Helmholtz problem rather than the biharmonic problem. Addition of this type of extra variability, no matter how small, shows that the paradox evaporates, i.e. the existence of the paradox is unstable. Villas Boas [@mario], who has considered the problem from a three-dimensional perspective, also demonstrates that the paradox evaporates when extra [*spatial*]{} variability is incorporated.

In subsequent analysis we will use the Laplace-Helmholtz model as a approximation for the time-[*independent*]{} 2D case, but where the modification is regarded as a simple way as approximating the combined effect of the non-linearities via a linear term. The question is as to whether we can do this in a sensible and self-consistent manner. We shall look at this in a variety of ways, starting with a rather [*ad hoc*]{} approach.

Reorganizing the low Reynolds number\
equations
=====================================

The production of a global solution to the equations requires a slightly unusual approach. We take the time-independent equation $$\nabla^4 \psi = R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)} \label{oldfour}$$ as our starting point and note some obvious facts. First, if we have a $\psi_0$ satisfying Laplace’s equation, then it satisfies equation Eq. (\[oldfour\]) identically. Second, in the high Reynolds number limit, although this is a singular limit, we note that the Jacobian term should then vanish identically, so that any solution of $$\delsq \psi = f(\psi)$$ satisfies Eq. (\[oldfour\]) when $R_e \rightarrow \infty$. This will be true of course in the linear case when $f(\psi) \propto \psi$ and this last equation is of Helmholtz type. Now let $\lambda$ be a (possibly complex) number of order 1 and $\beta$ a parameter to be determined. In particular, if we consider solutions, $\psi_1$, of a parametrized Helmholtz equation in the form $$\delsq \psi_1 = \lambda R_e^{2\beta} \psi_1$$ then the Jacobian will vanish and we note that $$\nabla^4 \psi_1 - R_e \frac{\partial(\psi_1, \delsq \psi_1)}{\partial(x, y)} \equiv \lambda^2 R_e^{4\beta} \psi_1$$ Next suppose that we consider a more general stream function $\psi$ as $$\psi = \psi_0 + \psi_1$$ Then a trivial calculation tells us that $$\nabla^4 \psi - R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)} = \lambda^2 R_e^{4\beta} \psi_1 - \lambda R_e^{2\beta+1} \frac{\partial(\psi_0, \psi_1)}{\partial(x, y)}   \label{expandedns}$$ Our approach is therefore to combine solutions for potential flow, $\psi_0$, with a solution $\psi_1$ of the Helmholtz equation Eq. (23) that also satisfies the fourth order equation given by Eq. (24) without the identically zero non-linear term: $$\nabla^4 \psi_1 = \lambda^2 R_e^{4\beta} \psi_1$$ When $R_e \rightarrow 0$ this PDE of course approaches the biharmonic equation provided $\beta >0$. Then Eq. (\[expandedns\]) offers the possibility that we can solve a PDE of the form: $$\nabla^4 \psi - R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)} = O(R_e^k) f(x,y)$$ where $k$ and $f$ can be calculated and analyzed. The introduction of the term of order $R_e^{4\beta}$ in the right side of Eq. (\[expandedns\]) may seem like an artificial device, but given that (a) we are seeking solutions of the problem as $R_e \rightarrow 0$; (b) we are not modifying higher derivatives in the equations; (c) we have a valid perturbation equation in Eq. (28), we consider that it is valid to proceed with this approach.

Other justifications
--------------------

The author recognizes that the argument in the previous sub-section is somewhat [*ad hoc*]{}, and initially gives us no idea how to justify the choice of $\lambda$ or $\beta$, or, as we shall introduce, the relevant composite parameter $\epsilon$. However, there is another other route to justifying this approach. Let $\omega = \nabla^2 \psi$. Then the full 2D time-independent incompressible Navier–Stokes equation is precisely $$\nabla^2  \omega- R_e \frac{\partial(\psi,  \omega)}{\partial(x,y)}=0\ .$$ We already know from the existence of the Stokes paradox that simple linearization by setting $R_e=0$ is hopeless. We might therefore consider starting from other linearizations, for example, the Helmholtz linearization $$\nabla^2 \omega - \epsilon^2 \omega=0\ ,$$ where $\epsilon$ has to be determined either from experiment or from deeper theoretical considerations. The linear starting point for analysis is then not the biharmonic system but the Helmholtz-Laplace equation: $$\nabla^4 \psi - \epsilon^2 \nabla^2\psi=0\ .$$ In principle then we can consider justifying the choice of $\epsilon$ by matching $ \epsilon^2 \omega$ with $R_e \frac{\partial(\psi, \omega)}{\partial(x,y)}$ as closely as possible under a suitable norm. We will prove in the next section that [*for all $\epsilon>0$ the Helmholtz-Laplace equation admits a solution satisfying both the boundary conditions on the circle and the velocity condition at infinity*]{}, obtained by considering elementary solutions that are linear combinations of solutions of the Laplace equation and solutions of the Helmholtz equation.

Note that this approach and our [*ad hoc*]{} scheme are linked by the relationship $$\epsilon^2 = \lambda R_e^{2\beta}$$ The determination of $\lambda$ and $\beta$ are interesting challenges. We shall see in Section 8 that there are good theoretical grounds for setting $\beta=1$. For now we will leave both parameters general.

The case of asymptotically uniform flow past a circle
=====================================================

We now consider the well-trodden route to the analysis of a stream function associated with a flow that is uniform at infinity and satisfies appropriate boundary conditions on $r=1$. Using polar coordinates, we therefore want $$\psi(r,\theta) \sim r \sin\theta  \ \ {\rm as}\ \  r \rightarrow \infty$$ and $$\psi = 0 = \frac{\partial \psi}{\partial r}\ \  {\rm on}\ \  r=1 \label{circlebc}$$ We build the solution for $\psi_0$ and $\psi_1$ as follows, under the assumption that $\psi_0$ satisfies the Laplace equation and $\psi_1$ the Helmholtz conditions. Any sum of the two will satisfy the Laplace-Helmholtz equation. Given that the angular behaviour at infinity is $\propto \sin\theta$ we make the standard assumption and seek solutions $$\psi_i =  g_i(r) \sin\theta\ \ {\rm for}\ \ i=0,1$$ The potential flow part (as usual) will be taken to be $$g_0 = 1 + B/r$$ for some constant $B$. Our analysis will differ from Stokes’ classic (Stokes, 1851) treatment in that $g_1$ does not satisfy the separated biharmonic equation. Instead we use the Helmholtz equation, which in separated form is just $$g_1''(r) +\frac{1}{r}g_1'(r)-\left(\epsilon^2+\frac{1}{r^2}\right)
   g_1(r)=0$$ where $\epsilon = \sqrt{\lambda}R_e^{\beta}$. The general solution to this ODE is given by $$g_1(r) = \alpha K_1(\epsilon r) + \beta I_1(\epsilon r)$$ where $I_1, K_1$ are standard “modified” Bessel functions. It is now quite clear that we can construct $g_1$ so as the preserve the boundary condition at infinity by setting $\beta = 0$. The $K_1$ function decays exponentially if $\epsilon>0$. We are left with two arbitrary constants that can be determined by satisfying the boundary conditions on the circle given by Eq. (\[circlebc\]). This is a matter of elementary algebra using some standard Bessel function identities. The final result for the total stream function can be simplified to: $$\psi = \left[r
-\biggl(1+
\frac{2
   K_1(\epsilon )}{\epsilon  K_0(\epsilon
   )}
\biggr) 
\frac{1}{r}
+\frac{2 K_1(r \epsilon )}{\epsilon 
 K_0(\epsilon )}
 \right] \sin (\theta ) \label{streamfn}$$ We note some interesting facts about this expression. First, it has the [*right*]{} behaviour as $r \rightarrow \infty$ provided $\epsilon >0$, as the Bessel function of $r$ decays exponentially. Second, if we keep $r$ fixed and finite and let $\epsilon \rightarrow 0_+$ we obtain $$\psi \sim \frac{\sin\theta}{\log(2/\epsilon) - \gamma}  \left(r \log r - \frac{1}{2}r + \frac{1}{2r}\right) \label{nearcyl}$$ and we recover a multiple of Stokes’ (1851) solution satisfying the boundary conditions on the sphere (but not at infinity), and, furthermore, the multiple is now reminiscent of that arising from the methods of matched asymptotic expansions. The claim therefore is that it is Eq. (\[streamfn\]) that essentially resolves the paradox, as all the boundary conditions are satisfied, but the low Reynolds number limit for finite $r$ does not. The third observation is that in the neighbourhood of the surface of the sphere, $$\psi \sim \sin\theta (r-1)^2 \frac{\epsilon K_1(\epsilon)}{K_0(\epsilon)}$$ and that for small $\epsilon$ this is given by $$\frac{(r-1)^2 \sin (\theta
   )}{-\log (\epsilon )+\log (2)-\gamma
   }+O\left(\epsilon ^2\right)$$ We also note that the limit as $\epsilon \rightarrow \infty$ is just the potential flow limit: $$\psi \sim (r - \frac{1}{r}) \sin \theta$$ [*Mathematica*]{} code for the evaluation of the stream function and velocity field is given in the Appendix.

Choice of $\lambda$ and $\beta$
-------------------------------

At this stage we have no knowledge about how to fix the parameter $\epsilon =  \sqrt{\lambda}R_e^{\beta}$, and we now turn our attention to this issue and a more detailed analysis of Eq. (\[expandedns\]). We note first that our working assumption of a power dependence of $\epsilon$ on $R_e$ has not actually been necessary. All we need is that $\epsilon \rightarrow 0_+$ as $R_e \rightarrow 0$. We also note that the choice $\beta=1$ and $\lambda = 1/(4 e)$ recovers the result (Van Dyke, 1964) that the limiting stream function for fixed $r$ and small $R_e$ is, from Eq. (\[nearcyl\]) $$\psi \sim \frac{\sin\theta}{\log(4/R_e) + \frac{1}{2} - \gamma}  \left(r \log r - \frac{1}{2}r + \frac{1}{2r}\right)$$ so that the low $R_e$ limit matches exactly the first term of the MAE result.[^5] Clearly other choices of $\lambda$ can be considered, as can other powers of $R_e$ or even a more general $\epsilon$ still. However, in Section 8 we will establish a theoretical basis for estimating these parameters, and argue that $\beta=1$.

Analysis of the remainder
=========================

We must now analyze the right side of Eq. (\[expandedns\])we call it $\Gamma$and we do so without any assumption as the the form of $\epsilon$. We have, changing variables in the Jacobian to polar coordinates: $$\Gamma = \epsilon^4 \psi^1 - R_e \frac{\epsilon^2}{r}\left[ \frac{\partial \psi_0}{\partial r}\frac{\partial \psi_1}{\partial \theta} -  \frac{\partial \psi_1}{\partial r}\frac{\partial \psi_0}{\partial \theta}  \right]$$ Recall that we set $\psi_i = g_i(r) \sin\theta$. Having imposed the boundary conditions, we have $$g_0 = r - \frac{1}{r}\left(1 + \frac{2 K_1(\epsilon)}{\epsilon K_0(\epsilon)}  \right), \ \ 
g_1 = \frac{2 K_1(\epsilon r)}{\epsilon K_0(\epsilon)}$$ It follows that $$\Gamma = \epsilon^4 \sin\theta g_1 - R_e \epsilon^2 \sin\theta \cos\theta J(\epsilon, r) = \Gamma_{HA} +\Gamma_{NL}$$ where $J$ is the reduced Jacobian: $$J = \frac{1}{r}\left(g_1 \frac{\partial g_0}{\partial r} - g_0 \frac{\partial g_1}{\partial r}\right)$$ and $\Gamma_{HA} =  \epsilon^4 \sin\theta g_1 $ denotes the ‘Helmholtz artifact’ introduced by our approach, and $\Gamma_{NL}$ denotes the non-linear Jacobian term. Some calculation with Bessel function identities leads to $$J(\epsilon, r) = \frac{2 \left(K_0(\epsilon ) K_2(r \epsilon )-\frac{1}{r^2}K_0(r \epsilon ) K_2(\epsilon)\right)}{K_0(\epsilon )^2}$$ from which it is manifest that $J(\epsilon,1) = 0$ and the non-linear terms vanish on the cylinder boundary irrespective of the choice of $\epsilon$. Furthermore, the asymptotic behaviour of the Bessel functions tells us that provided $\epsilon>0$, $$J(\epsilon, r) = O(\frac{1}{\sqrt{r}} \exp(-\epsilon r))$$ as $r \rightarrow \infty$. We deduce that $\Gamma_{NL}$ is a bounded function for all $r$. What is its order as a function of $R_e$? A straightforward estimate may be given by looking in a neighbourhood of the cylinder. It is straightforward to establish that $$\frac{\partial J}{\partial r}|_{r=1} = \left(\frac{2 K_1(\epsilon)}{K_0(\epsilon)}  \right)^2 \sim \frac{4}{\epsilon ^2 (-\log (\epsilon )+\log (2)-\gamma )^2}$$ where the latter expression applies as $\epsilon \rightarrow 0$. We deduce that in the [*immediate*]{} neighbourhood of the cylinder $$\Gamma_{NL} = O\left(\frac{R_e}{(\log\epsilon)^2}\right)$$ which is $o(R_e)$ provided only that $\epsilon \rightarrow 0$ as $R_e \rightarrow 0$. Although $J$ grows to a maximum as $r$ is increased away from unity, before reaching a maximum and then decaying for large $r$, some numerical experiments with [*Mathematica*]{} confirm that the maximum of $\epsilon^2 J$ slowly decreases as $\epsilon$ decreases to zero. So we can assert that $\Gamma_{NL}$ is well behaved and is $o(R_e)$. When we consider $\Gamma_{HA}$, it is easy to see that $g_1$ has a maximum on $r=1$, and that $$\Gamma_{HA}|_{r=1} = \epsilon^4 \sin\theta g_1(1) = O(\epsilon^2/\log\epsilon)$$ as $\epsilon \rightarrow 0$. So this term also behaves. If we desire that the Helmholtz artifact tends to zero faster than the non-linear term (which is desirable for the credibility of our approach) it is then natural to specify $\epsilon$ in the form $\sqrt{\lambda}R_e^{\beta}$ and to demand that $\beta \ge 1$. In particular, the choice $\beta=1$ arranges that $$\Gamma_{HA} =  O(R_e^2/\log R_e)\ ,\ \ \Gamma_{NL} = O(R_e/\log^2{R_e})$$ and we have established that $\Gamma = o(R_e)$. It also decays exponentially at infinity.

The drag on the cylinder
========================

The calculation of the drag can be done by integrating the pointwise force on the cylinder over its surface. The pointwise force has two components. One involves the local rate of strain in the fluid, and the other is the pressure force. The first requires a purely local calculation, but the second requires an integration of the pressure equation from infinity to the cylinder. A question is how this second part can be carried out [*without*]{} any [*global*]{} representation of the flow field. For the [*specific*]{} case of low Reynolds number calculations with a certain symmetry there are ways around the problem that we shall exploit presently. To make these matters clear, we shall summarize a first-principles calculation of the force using the results for cylindrical polar coordinates for the rate of strain as given by Appendix A of Acheson (1990). The fluid velocity is given (in our units) by $$\underline{u} = -\frac{U}{r}\frac{\partial \psi}{\partial \theta} \underline{e}_r + U \frac{\partial \psi}{\partial r} \underline{e}_{\theta}$$ On the cylinder the rate-of-strain tensor has components $e_{rr} = 0 = e_{\theta \theta}$, and $$2e_{r\theta} = U \frac{\partial^2 \psi}{\partial r^2}$$ The stress tensor $T_{ij} = -p \delta_{ij} + 2\mu e_{ij}$, so the force on the cylinder boundary is $$\underline{t} = T.\underline{e}_r =  -p \underline{e}_r  + \mu U \frac{\partial^2 \psi}{\partial r^2} \underline{e}_{\theta}$$ We project this into the $x$-direction and integrate over the circle to get the following formula for the drag $D_0$ (force per unit length on the cylinder): $$D_0 = - \int_0^{2\pi}\left[p \cos\theta + \mu U \frac{\partial^2 \psi}{\partial r^2} \sin\theta \right]d\theta$$ For low $R_e$ we estimate the pressure using the momentum equation in the form $$\del p = -\mu \del \wedge (\del \wedge \uvec )$$ and a short calculation gives $$\del p = \mu U \left\{-\frac{1}{r} \frac{\partial\ }{\partial \theta}(\delsq \psi)\ervec + \frac{\partial\ }{\partial r}(\delsq \psi)\ethvec \right\}$$ We need to integrate this in from infinity to a general point on the circle $r=1$. In general this is awkward without a global form of $\psi$. But if $\psi$ has the form $\psi = g(r) \sin\theta$ then we have $$\delsq \psi = \sin\theta {\cal L}[g(r)]\  \ {\rm where}\ \ {\cal L} [g(r)] \equiv  \frac{1}{r}\frac{\partial\ }{\partial r}\left(r\frac{\partial g}{\partial r} \right)- 
\frac{g}{r^2}$$ and the pressure equation becomes: $$\del p = \mu U \left\{-\ervec \frac{1}{r} \cos\theta {\cal L}[g(r)] + \ethvec \sin\theta \frac{\partial\ }{\partial r}{\cal L}[g(r)]\right\}$$ We do the integration from infinity along $\theta = \pi/2$ and then work around the circle $r=1$. Carrying this out we obtain $$p(1,\theta) = p_{\infty} - \mu U \cos \theta \frac{\partial\ }{\partial r} {\cal L}[g(r)]|_{r=1}$$ We now evaluate the total integral for the force per unit length to obtain $$D_0 = \mu U \pi \left[\frac{\partial\ }{\partial r}{\cal L}[g(r)]-\frac{\partial^2 \!g(r)}{\partial r^2} \right]$$ The drag coefficient for such a 2D problem is defined following Tritton (1959) as $C_D = |D_0/(\rho U^2)|$ and is now easily seen to be given by $$C_D = \frac{\pi}{R_e} \left|\frac{\partial\ }{\partial r}{\cal L}[g(r)]-
\frac{\partial^2\!g(r)}{\partial r^2} \right|$$ evaluated on $r=1$. Given the vanishing of $g(1)$ and $g'(1)$ this simplifies further to $$C_D = \frac{\pi}{R_e} |g^{(3)}(1) |$$ In the case given by the MAE approach the stream function is given in the neighbourhood of the cylinder by Eq. (1) and then $g$ and $C_D$ are given by $$g(r) = C\left(r \log r -\frac{r}{2} + \frac{1}{2r}\right),\ \ C_D = \frac{4 \pi C}{R_e}$$ In our new model $g$ as is given by the radial part of Eq. (\[streamfn\]), and this time, using some Bessel identities, we obtain $$C_D = \frac{2 \pi}{R_e} \frac{\epsilon^2 K_2(\epsilon)}{K_0(\epsilon)}$$

![Comparison with Tritton (1959) data](tritton){width="3.0in"}

We recall now the original assumption that $\epsilon = \sqrt{\lambda}R_e^{\beta}$. If we fix $\beta = 1$ but allow $\lambda$ to vary the drag coefficient of our model is then $$C_D = 2 \pi \lambda R_e \frac{K_2(\sqrt{\lambda}R_e)}{K_0(\sqrt{\lambda}R_e)}$$ We can plot the drag coefficient from our model, assuming that $\beta = 1$, with various choices for $\lambda$, and compare the results with those of Tritton (1959), which used an $R_e$ based on the diameter. Tritton’s data has been converted by halving [*his*]{} Reynolds number. Also, in making a comparison with Fig. 8.5 of Van Dyke (1964), it appears that Van Dyke has plotted $C_D/(4 \pi)$, as otherwise it is not possible to reconcile that plot with Tritton’s data. In Figure 1 we show the data for the first three fibres used by Tritton (those that consider the lowest $R_e$) and

-   our model with $\lambda$ picked to match the MAE result, $\lambda = 1/(4 e) \sim 0.09197$;

-   our model with a least squares best fit $\lambda \sim 0.04452$;

-   the one term MAE result based on $ C = \Delta_1(R_e) = (\log(4/R_e) - \gamma + 1/2)^{-1}$;

-   the two term MAE result based on $ C = \Delta_1(R_e)- 0.8669 \Delta_1^3(R_e)$.

The choice of constants in the MAE approach is itself somewhat arbitrary. That made by Kaplun (1957), giving Eq. (\[nearcyl\]), is little more than convenience. These results give support to the new model and the agreement with experiment we get by taking $\lambda$ about one half that implied by the MAE approach is rather tantalizing.

Iteration: a Stokes-Whitehead anomaly?
======================================

Another use for the type of low Reynolds numbers solution developed here is to provide a basis for [*iteration*]{} of the solution in powers of $R_e$. The presence of the paradox obstructing a global representation in the standard approach makes this impossible. Furthermore, we know from a corresponding analysis of the 3D spherical problem that even when the base solution makes sense, even the first iteration may fail: the [*Stokes-Whitehead*]{} paradox emerges. So it is of considerable interest to investigate the iteration of the system. In doing so we will keep $\epsilon$ general initially.

We summarize the solution developed thus far. With $\psi$ as given we have, combining equations (\[expandedns\],45) and our definitions , $$\nabla^4 \psi - R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)} = \Gamma = \Gamma_{HA}+\Gamma_{NL}\ ,$$ where $$\Gamma_{HA} = \epsilon^4 \sin\theta g_1   =  \epsilon^4 \sin\theta \frac{2 K_1(\epsilon r)}{\epsilon K_0(\epsilon)} =\gamma_{HA}(\epsilon,r)  \sin\theta$$ and $$\begin{split}
\Gamma_{NL} &=  - R_e \epsilon^2 \sin\theta \cos\theta J(\epsilon, r)\\ &=- R_e \epsilon^2 \sin\theta \cos\theta  \frac{2 \left(K_0(\epsilon ) K_2(r \epsilon )-\frac{1}{r^2}K_0(r \epsilon ) K_2(\epsilon)\right)}{K_0(\epsilon )^2}\\
&= \gamma_{NL}(\epsilon, r)\sin\theta \cos\theta 
\end{split}$$ where these equations also serve to define the radial terms $\gamma_{XY}$. From now on we will work in terms of these radial functions, so it is useful just to write down their explicit forms: $$\gamma_{HA}(\epsilon,r) = 2 \epsilon^3 \frac{K_1(\epsilon r)}{K_0(\epsilon)}$$ $$\gamma_{NL}(\epsilon,r) = -2 R_e \epsilon^2 \frac{ \left(K_0(\epsilon ) K_2(r \epsilon )-\frac{1}{r^2}K_0(r \epsilon ) K_2(\epsilon)\right)}{K_0(\epsilon )^2}$$

Development of the iteration
----------------------------

We now attempt to refine the solution by writing the total solution as $$\Psi = \psi + \psi_2$$ Ideally we would like to arrange that $$0 = \nabla^4 \Psi - R_e \frac{\partial(\Psi, \delsq \Psi)}{\partial(x, y)}$$ The right side of this is given by the expansion $$\nabla^4 \psi - R_e \frac{\partial(\psi, \delsq \psi)}{\partial(x, y)}+\nabla^4 \psi_2 - R_e \frac{\partial(\psi_2, \delsq \psi)}{\partial(x, y)}- R_e \frac{\partial(\psi, \delsq \psi_2)}{\partial(x, y)}
- R_e \frac{\partial(\psi_2, \delsq \psi_2)}{\partial(x, y)}$$ and we write the remaining parts of the Navier-Stokes equation as $$\begin{split}
0 = & \nabla^4 \psi_2  +\gamma_{NL}(\epsilon, r)\sin\theta \cos\theta +\gamma_{HA}(\epsilon,r)  \sin\theta \\
& - R_e \frac{\partial(\psi_2, \delsq \psi)}{\partial(x, y)}- R_e \frac{\partial(\psi, \delsq \psi_2)}{\partial(x, y)}
- R_e \frac{\partial(\psi_2, \delsq \psi_2)}{\partial(x, y)}
\end{split}$$ We do not know how to solve this full non-linear system, and propose instead to treat a linearized form. The question now is whether to work with the raw system $$0 =  \nabla^4 \psi_2  +\gamma_{NL}(\epsilon, r)\sin\theta \cos\theta +\gamma_{HA}(\epsilon,r)  \sin\theta \\$$ or to propose a Helmholtz-damped variation along the lines of our original approach. We also need to understand which non-linear term to treat first. [*Provided*]{} $\epsilon$ is of the form $\epsilon = \sqrt{\lambda} R_e$, for example, we know that $\Gamma_{NL}$ is the lowest order in $R_e$: $\Gamma_{NL} = O(R_e^2/\log^2 R_e)$, so we will give priority to killing this lowest order correction.

Naive analysis of the raw forms
-------------------------------

To develop the solution we need to identify the Greens’ s function for the biharmonic operator. In fact we do not need the full form as the angular structure of the right side is simple - we need only look for a pair of appropriate radial Green’s functions satisfying appropriate boundary conditions on the cylinder and at infinity. To this end we write $$\psi_2(r, \theta) = \psi_{3}(r) \sin\theta \cos\theta + \psi_4(r) \sin \theta$$ and require that $$\nabla^4 (\psi_{3}(r) \sin\theta \cos\theta) = - \gamma_{NL}(\epsilon, r)\sin\theta \cos\theta$$ and $$\nabla^4 (\psi_{4}(r) \sin\theta) = - \gamma_{HA}(\epsilon, r)\sin\theta$$ We also require that $\psi_i$ vanishes as $r \rightarrow \infty$ and both $\psi_i$ and its first derivative vanish on $r=1$. To proceed further we recall the elementary form of the Laplacian in 2D cylindrical polar coordinates: $$\nabla^2 f(r, \theta) = \frac{\partial^2 f}{\partial r^2} + \frac{1}{r} \frac{\partial f}{\partial r} + \frac{1}{r^2}  \frac{\partial^2 f}{\partial \theta^2}$$ So if $f$ is of the form $f = h(r) \sin(k \theta)$ then $$\nabla^2 f(r, \theta) = \sin(k \theta) \biggl(\frac{d^2 h}{d r^2} + \frac{1}{r} \frac{dh}{dr} - \frac{k^2}{r^2} h\biggr)$$ Define the $k$’the radial modal Laplacian operator as $${\cal L}_k = \frac{d^2 \ }{d r^2} + \frac{1}{r} \frac{d\ }{dr} - \frac{k^2}{r^2}$$ We now have a pair of ODEs in the form $${\cal L}_2^2 \psi_3(r) =  - \gamma_{NL}(\epsilon, r)$$ $${\cal L}_1^2 \psi_4(r) =  - \gamma_{HA}(\epsilon, r)$$ In a simpler problem we would now construct a pair of Green’s functions $G_i(r,x)$, $i=1,2$, with the properties that $${\cal L}_i^2 G_i(r,x) = \delta(r-x)$$ subject to, [*if it is possible*]{}, $$\lim_{r \rightarrow \infty}G_i(r,x) = 0$$ $$G_i(1,x) = 0$$ $$\frac{\partial\ }{\partial r}G_i(r,x)|_{r=1} = 0$$ However, there are significant complications. The Helmholtz term we have thought of as an modelling artifact has first to be understood as an optimal linear approximation to the full non-linear system. The first steps on this are considered in Section 8. We will eventually have to manage the fact that no $G_1$ function exists satisfying these boundary conditions. However, it is now consistent to proceed to construct $G_2$, as this (a) exists, and (b) is needed for the management of the full non-linear system and its solution.

Integration of the inertia term
-------------------------------

We will now focus on the lowest order correction arising from the inertia term, for which $G_2$ is the required radial Green’s function. The building blocks for this are two different solutions of the ODE $$\biggl(\frac{d^2 \ }{d r^2} + \frac{1}{r} \frac{d\ }{dr} - \frac{2^2}{r^2}\biggr)^2 \phi(r) = 0\ .$$ The solution to the homogeneous problem is of the form $$\phi(r) = a r^4 + b r^2 + c + \frac{d}{r^2}$$ To try to build a Green’s function, we would write down a pair of solutions: $$\begin{split}
\phi_L(r) &= a_L r^4 + b_L r^2 + c_L + \frac{d_L}{r^2}\ ,\ \ 1 \leq r \leq x\ ,\\
\phi_R(r) &= a_R r^4 + b_R r^2 + c_R + \frac{d_R}{r^2}\ ,\ \ x \leq r < \infty\ .
\end{split}$$ The condition at infinity requires that we set $a_R = 0 = b_R$, so that $$\phi_R(r) =  c_R + \frac{d_R}{r^2} ,\ \ x \leq r < \infty\ .$$ The boundary conditions on the cylinder require that $$\begin{split}
0 &= a_L+b_L+c_L+d_L \\
0 &= 4a_L + 2 b_L-2d_L
\end{split}$$ So the inner function is given by $$\phi_L(r)= \frac{2 a_L+b_L}{r^2}+a_L r^4+b_L r^2-2 b_L-3 a_L$$ We now have four unknowns that must be obtained by imposing junction conditions to obtain the delta-function. The junction conditions are: $$\begin{split}
\phi_R(x) &= \phi_L(x)\\
\phi'_R(x) &= \phi'_L(x)\\
\phi''_R(x) &= \phi''_L(x)\\
\phi_R^{(3)}(x) &= \phi_L^{(3)}(x)+1
\end{split}$$ These equations were solved and simplified using [*Mathematica*]{} V7 and the solutions are $$\begin{split}
a_L &= --\frac{1}{48 x}\\
b_L &= \frac{1}{16}x\\
c_R &= \frac{\left(x^2-1\right)^2}{16x}\\
d_R &= -\frac{x^6-3 x^2+2}{48 x}
\end{split}$$ Finally the two parts of the Green’s function are (introducing the $x$-dependence explicitly) $$\begin{split}
\phi_L(r,x) &=-\frac{r^4}{48 x}+\frac{r^2 x}{16}+\frac{\frac{x}{16}-\frac{1}{24x}}{r^2}-\frac{x}{8}+\frac{1} {16 x} \ ,\ \ 1 \leq r \leq x\ ,\\
\phi_R(r,x) &= -\frac{(x^6-3 x^2+2)}{48 r^2x}+\frac{\left(x^2-1\right)^2}{16x}\ ,\ \ x \leq r < \infty\ .
\end{split}$$ The solution for $\psi_3$ is then given formally by $$\psi_3(r) = -\int_1^r \gamma_{NL}(x) \phi_R(r,x)dx -\int_r^{\infty} \gamma_{NL}(x) \phi_L(r,x)dx$$ This stream function is $O(1)$ as $r \rightarrow \infty$ so we claim that a secondary paradox of Stokes-Whitehead type has been avoided. To proceed further we need to better understand the full non-linear system.

The non-linear equations and constraints on $\epsilon$
======================================================

The task now at hand is to develop a proper theoretical basis for the estimation of the parameter $\epsilon$. To this end we must write down the full equations. Converting to polar coordinates, the Navier–Stokes equations under consideration are $$\nabla^4 \psi = \frac{R_e}{r} \frac{\partial(\psi, \delsq \psi)}{\partial(r, \theta)}$$ Let us assume that $R_e$ is small enough that the flow remains symmetric about the horizontal axis. We may then write a Fourier decomposition $$\psi = \sum_{m=1}^{\infty}\sin(m\theta) \phi_m(r)$$ Then with ${\cal L}_k$ as above, $$\nabla^2 \psi = \sum_{m=1}^{\infty}\sin(m\theta) {\cal L}_m \phi_m(r)$$ $$\nabla^4 \psi = \sum_{m=1}^{\infty}\sin(m\theta) {\cal L}_m^2 \phi_m(r)$$ $$\begin{split}
&\sum_{m=1}^{\infty}\sin(m\theta) {\cal L}_m^2 \phi_m(r) \\ &=  \frac{R_e}{r} \sum_{l=1}^{\infty} \sum_{n=1}^{\infty} \biggl[\sin(l\theta)n \cos(n\theta)\frac{\partial \phi_l}{\partial r}{\cal L}_n \phi_n - l \cos(l\theta) \sin(n\theta) \phi_l \frac{\partial \ }{\partial r}{\cal L}_n \phi_n  \biggr]\\
\\ &=  \frac{R_e}{2 r} \sum_{l=1}^{\infty} \sum_{n=1}^{\infty} \biggl[n[\sin((l-n)\theta)+\sin((l+n)\theta)]\frac{\partial \phi_l}{\partial r}{\cal L}_n \phi_n \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ - l [\sin((l+n)\theta)-\sin((l-n)\theta)]\phi_l \frac{\partial \ }{\partial r}{\cal L}_n \phi_n  \biggr]
\end{split}$$ Doing the Fourier analysis gives us $$\begin{split}
&{\cal L}_m^2 \phi_m(r) \\ =  \frac{R_e}{2r}  \biggl[& (1-\delta_{m,1})\sum_{l=1}^{m-1}\biggl((m-l)\frac{\partial\phi_l}{\partial r}{\cal L}_{m-l}\phi_{m-l} - l \phi_l \frac{\partial \ }{\partial r}{\cal L}_{m-l}\phi_{m-l}\biggr)\\
&+\sum_{l=1}^{\infty}\biggl(l \frac{\partial \phi_{l+m}}{\partial r}{\cal L}_l \phi_l+(l+m)\phi_{l+m}\frac{\partial\ }{\partial r}{\cal L}_{l}\phi_{l}\biggr) \\
&-\sum_{l=1}^{\infty}\biggl((l+m) \frac{\partial \phi_l}{\partial r}{\cal L}_{l+m} \phi_{l+m}+l\phi_l\frac{\partial\ }{\partial r}{\cal L}_{l+m}\phi_{l+m}   \biggr)\biggr]
\end{split}$$ Thus far no approximations have been made, other than to assume the flow remains symmetric about the $x$-axis. Now we approximate the model by considering only the contributions of the terms with $m=1,2$. We then have the simpler, but coupled and non-linear system $${\cal L}_2^2 \phi_2 = \frac{R_e}{2r}\biggl(\frac{\partial \phi_1}{\partial r}{\cal L}_1\phi_1 - \phi_1 \frac{\partial\ }{\partial r}{\cal L}_1 \phi_1  \biggr) \label{phitwonl}$$ $${\cal L}_1^2 \phi_1 = \frac{R_e}{2r}\biggl(\frac{\partial \phi_2}{\partial r}{\cal L}_1\phi_1 - \phi_1 \frac{\partial\ }{\partial r}{\cal L}_2 \phi_2 -2\frac{\partial \phi_1}{\partial r}{\cal L}_2\phi_2 +2 \phi_2 \frac{\partial\ }{\partial r}{\cal L}_1 \phi_1 \biggr)\label{phionenl}$$ With our simple form of $\phi_1$ that is a solution of the Laplace-Helmholtz condition $${\cal L}_1^2 \phi_1 = \epsilon^2 {\cal L}_1 \phi_1\label{philin}$$ we have already determined the solution of the Eq. (\[phitwonl\])– $\phi_2$ is just $\psi_3$ as already given above. In order to generate maximal self-consistency of the linearization we must choose $\epsilon$ in such a way as to minimize the mismatch between the right sides of Eq. (\[phionenl\]) and Eq. (\[philin\])[^6]

Let us look at the overall scaling behaviour in terms of the Reynolds number. We let $\phi_2 = R_e \tilde{\phi}_2$, so that $${\cal L}_2^2 \tilde{\phi}_2 = \frac{1}{2r}\biggl(\frac{\partial \phi_1}{\partial r}{\cal L}_1\phi_1 - \phi_1 \frac{\partial\ }{\partial r}{\cal L}_1 \phi_1  \biggr)$$ The required matching is then of the form $$\epsilon^2 {\cal L}_1 \phi_1 \sim {\cal L}_1^2 \phi_1 = \frac{R_e^2}{2r}\biggl(\frac{\partial \tilde{\phi}_2}{\partial r}{\cal L}_1\phi_1 - \phi_1 \frac{\partial\ }{\partial r}{\cal L}_2 \tilde{\phi}_2 -2\frac{\partial \phi_1}{\partial r}{\cal L}_2\tilde{\phi}_2 +2 \tilde{\phi}_2 \frac{\partial\ }{\partial r}{\cal L}_1 \phi_1 \biggr)$$ which strongly supports the scaling behaviour $$\epsilon^2 = \lambda R_e^2$$ as $R_e \rightarrow 0$. That is, we have a theoretical basis for setting $\beta=1$. A determination of $\lambda$ requires the introduction of a suitable rigorous criteria for minimizing the mismatch between the linear and non-linear forms, and this is under investigation.

Note that the [*full*]{} equations, but limited to the first two angular modes, can certainly be written [*without further approximation*]{} as $${\cal L}_1^2 \phi_1 = R_e^2 \lambda(r) {\cal L}_1 \phi_1$$ for some unknown function $\lambda(r)$, and our method can now be properly understood as that of working with some “average” value of $\lambda(r)$, and noting that the resulting solution is free of a paradox. Note also that writing the right side of this linearized system as a multiple of ${\cal L}_1 \phi_1$ is not as arbitrary as it might seem, for we know that the right side of the full non-linear system vanishes identically when ${\cal L}_1 \phi_1=0$. That is, this Laplace-Helmholtz model has a proper theoretical justification, rather than merely being the basis of a convenient interpolation between the boundary conditions on the cylinder and at infinity.

We could also consider generalizations where an improved [*ansatz*]{} for the form of $\lambda(r)$ is employed. The optimal average form and improved functional choices are under investigation. Only once this has been done would it make sense to consider further iteration. Some initial considerations suggest that as a function $\lambda(1)=0$ and that $\lambda(r)$ might be asymptotic to a constant independent of $R_e$ as $r \rightarrow \infty$, but further analysis is needed.

Summary
=======

We have constructed a [*global*]{} stream function satisfying the two-dimensional viscous incompressible steady equations in the limit $R_e \rightarrow 0$. The stream function satisfies the correct boundary conditions on a cylinder and infinity. The “perturbation” to the linearized equations introduced to accomplish this is of higher order in $R_e$ than the inertia terms, and this modification is now properly understood as a linear representation of the full non-linear theory. The results suggest improved agreement with experimental data over those obtained by the MAE approach. Further work is needed on this approach, in particular on comparisons with newer data sets and determination of the remaining free parameter $\lambda$ from theoretical considerations. Although the approach of this paper initially started with rather [*ad hoc*]{} considerations, the method developed here is founded on a deeper consideration based on the approximation of the non-linear equations with an optimal linear approximation based on equations of Laplace-Helmholtz type, rather than the biharmonic equation. The need to properly treat the non-linearities has also been illustrated by the observation that the presence of a paradox is unstable with respect to small changes in the system.

Model summary
-------------

The viscous and “paradox-free” stream function is given in non-iterated form by $$\psi = \left[r
-\biggl(1+
\frac{2
   K_1(\varepsilon )}{\varepsilon  K_0(\varepsilon
   )}
\biggr) 
\frac{1}{r}
+\frac{2 K_1(r \varepsilon )}{\varepsilon 
 K_0(\varepsilon )}
 \right] \sin (\theta )$$ and the parameter $\epsilon = \sqrt{\lambda}R_e$. The drag coefficient with Tritton’s conventions is: $$C_D = \frac{2 \pi}{R_e} \frac{\epsilon^2 K_2(\epsilon)}{K_0(\epsilon)}$$ The choice of $\lambda$ from theoretical considerations has yet to be determined, but experimental drag data suggests $\lambda \sim 0.04452$. The stream function satisfies the correct boundary conditions on the cylinder and at infinity for all $\epsilon>0$, and also reduces to that for potential flow as $\epsilon \rightarrow \infty$. [*Mathematica*]{} code for the stream function and velocity field are given in the Appendix.

Credit
------

I am grateful to Bin Zhou for his comments on the earlier (2006) version of this paper.

[10]{}

Acheson, D.J. 1990. [*Elementary Fluid Dynamics*]{}. Clarendon Press.

Kaplun, S. 1957 Low Reynolds number flow past a circular cylinder. [*J. Math. Mech.*]{} [**6**]{}, 595-603.

Lamb, H, 1911 On the uniform motion of a sphere through a viscous fluid. [*Philos. Mag.*]{} \[6\], [**21**]{}, 112-121.

Milne-Thomson, L.M. 1996 [*Theoretical Hydrodynamics*]{} Fifth Edition, Dover Publications.

Ockendon, J.R., Howison, S.D., Lacey, A. & Movchan, A. 2003 [*Applied Partial Differential Equations*]{} rev. ed. Oxford University Press.

Oseen, C.W. 1910, [*Ark. Math. Astronom. Fys*]{} [**6**]{}, No. 29.

Proudman, I., and Pearson, J.R.A. (1957). Expansions at small Reynolds numbers for the flow past a sphere and a circular cylinder. [*J. Fluid Mech.*]{} [**2**]{}, 237-262.

Stokes, G.G. 1851 On the effect of the internal friction of fluids on the motion of pendulums [*Trans. Camb. Phil. Soc.*]{} [**9**]{}, Part II, 8-106.

Tritton, D.J. (1959) Experiments on the flow past a circular cylinder at low Reynolds number. [*J. Fluid Mech.*]{} [**6**]{}, 547-567.

Van Dyke, M. 1964 [*Perturbation Methods in Fluid Mechanics*]{}, Academic Press.

Villas Boas, M.B., 2007. Stokes Paradox; A new approach (in Portuguese). Presentation at 8th Congreso Iberoamericano de Ingenieria Mecánica.

Appendix: [*Mathematica*]{} code for the flow field {#appendix-mathematica-code-for-the-flow-field .unnumbered}
===================================================

The following code may be useful. First the stream function:

    \[Psi][r_, \[Theta]_, \[Epsilon]_] := Sin[\[Theta]] (r - (1 + 
           2 BesselK[1, \[Epsilon]]/\[Epsilon]/BesselK[0, \[Epsilon]])/
         r + 2 BesselK[1, r \[Epsilon]]/\[Epsilon]/BesselK[0, \[Epsilon]])

    The streamlines are then easily visualized:
    ContourPlot[\[Psi][Sqrt[x^2 + y^2], ArcTan[x, y], 0.01], {x, -10 , 
      10}, {y, -5, 5}, AspectRatio -> 1/2, Contours -> 50, 
     RegionFunction -> Function[{x, y}, x^2 + y^2 >= 1], 
     Epilog -> Circle[{0, 0}, 1]]
     

![Streamlines](stream.pdf){width="3.5in"}

The Cartesian components of the velocity field are given by

     CartesianVelocity[x_, y_, U_, \[Epsilon]_] := 
     Module[{r = Sqrt[x^2 + y^2], \[Theta] = ArcTan[x, y], Ur, U\[Theta]},
       Ur = -((U*((-1 + r^2)*\[Epsilon]*BesselK[0, \[Epsilon]] - 
              2*BesselK[1, \[Epsilon]] + 2*r*BesselK[1, r*\[Epsilon]])*
            Cos[\[Theta]])/(r^2*\[Epsilon]*BesselK[0, \[Epsilon]]));
      U\[Theta] = 
       U*(1 + r^(-2) - (BesselK[0, 
              r*\[Epsilon]] - (2*
                BesselK[1, \[Epsilon]])/(r^2*\[Epsilon]) + 
             BesselK[2, r*\[Epsilon]])/BesselK[0, \[Epsilon]])*
        Sin[\[Theta]];
      Ur*{Cos[\[Theta]], Sin[\[Theta]]} + 
       U\[Theta]*{-Sin[\[Theta]], Cos[\[Theta]]}]
     

The vector flow field is shown in Figure 3:

     VectorPlot[
     If[x^2 + y^2 >= 1, 
      CartesianVelocity[x, y, 1, 0.1], {0, 0}], {x, -10, 10}, {y, -5, 5}, 
     RegionFunction -> Function[{x, y}, x^2 + y^2 >= 1], 
     AspectRatio -> 1/2 , VectorScale -> Small, 
     Epilog -> Circle[{0, 0}, 1]]
     

![Vector field of flow](vectors.pdf){width="3.5in"}

This illustrates the satisfaction of the boundary conditions on the cylinder.

[^1]: Corresponding author: Department of Mathematics King’s College, The Strand, London WC2R 2LS, England; E-mail: william.shaw@kcl.ac.uk

[^2]: This version April 2009, with a first consideration of iteration of the solution.

[^3]: Here $\delsq$ acting on [*vectors*]{} should be understood as the ordinary Laplacian acting on Cartesian components.

[^4]: it is also common to base $R_e$ on the diameter.

[^5]: for analyses basing $R_e$ on the diameter we have $\log(8/R_e)$ etc. in the denominator.

[^6]: This is of course one of many ways of proceeding, but our goal here is to provide a theoretical basis for the estimation of the so far free parameter $\epsilon$.
