Creeping flow past a stationary sphere

This notebook has been written in Mathematica by 

Mark J. McCready
Professor of Chemical Engineering
University of Notre Dame
Notre Dame IN 46556
USA

mjm@nd.edu

http://www.nd.edu/~mjm/

http://www.chemeprof.com


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  10/24/17

GITaylor movie:  https://www.youtube.com/watch?v=51-6QCJTAjU

This notebook shows how to solve creeping (inertialess) flow past a stationary sphere (Stokes's Problem)

Reference:  S. Middleman  (1999) An Introduction to Fluid Dynamics:Principles of Analysis and Design,
Wiley pp 166-171.

The flow situation of interest is shown in this figure.  Note that fluid motion past a stationary sphere is described with vr and vθ.  You would use vϕ if there were complete rotation.  Also think about what vr and vθ are as r-> infinity.  What trig function of θ will give behavior consistent with the orientation of the flow and the direction that θ is changing.  As a hint, at θ = 0, π, the flow is purely radial.  At θ = π/2, 3π/2 the flow is purely tangential.

Keywords

Creeping flow, Viscous flow theory, Stokes Flow, Low Reynolds number,Transport phenomena, Fluid Dynamics

Problem of interest

An issue that occurs in many situations of interest to chemical engineers is the contacting of solid particles with liquids or the separation of particles from gases and liquids.  In all of these situations, there is a need to either keep the particles suspended, or, to hasten their settling.  To understand this problem, we need to find an expression that gives the steady -state drag on a particle.  There is no simple relation that is valid for all Reynolds numbers.  However, if the particles is a sphere and the Reynolds number sufficiently small, a relation called "Stokes Law" is valid and it can be obtained through solution of the Navier-Stokes equations.

We will derive this relation below and in the process will hopefully learn something about the physics of "creeping" flows.  It is worthwhile to emphasize that even if a particle is not a sphere, the relation will give the qualitative prediction of the drag and should not be too wrong quantitatively.

Consider the flow past a stationary sphere where the Reynolds number is significantly less than unity.  Inertia can be neglected.  The far away velocity field is a constant straight flow.  To make the problem as simple as possible, we will use spherical coordinates with the ϕ axis oriented parallel to the flow.  Thus we have ϕ symmetry and no uϕ.  This means that the (Subscript[u, ϕ]) equation is not needed.  If we used a different coordinate system, or different orientation, the problem would be more complicated.

Learning Objectives

Physical Issues

1.  We are restricting the problem to the case where inertia forces are much weaker than viscous forces.  Thus we expect that the inertia terms of the Navier-Stokes equation should be much smaller than the viscous and pressure terms and that thus they can be neglected.  For this case the Reynolds number is very small.  If we make the equations dimensionless all terms are no larger in magnitude than about unity.  Thus the parameters that appear in these equations, and which can have values much different from on, determine which terms are needed for the solution.  In the nondimensional equations, the Reynolds number multiplies the inertia terms and these will consequently be neglected in the solution.  

2.  Since the Reynolds number is small, the fluid goes only where specifically pushed and it will stop if the forcing is stopped.  Thus the geometry of the flow field is determined by the boundaries.  

3.  Because viscous forces dominate the flow field, the fluid can never accelerate above the free stream value even if an obstacle causes the fluid to be squeezed.  Thus the velocity in the region of the sphere just slows down and then returns to the free stream value.  

4.  Both normal stresses and tangential stresses contribute to the drag on the sphere.  These can be termed form drag and skin drag.     

5.  Consistent with the fluid not accelerating, the pressure never increases above the free stream value.  The fluid has no inertia that would cause a pressure increase as the fluid slows down.

6.  The velocity decays slowly (as  (1/r)) and thus the disturbance is felt very far away from the sphere.  This makes it difficult to do a real experiment, in a reasonable size container, that allows that sphere to fall at a speed specified by the drag that is predicted from the analysis here.  The very high Reynolds number case decays much faster.

Mathematical Issues

1.  Linear and nonlinear partial differential equations:  In general nobody can make much progress solving nonlinear partial differential equations analytically.  Thus, if we are to solve problems involving the Navier Stokes equations, there will have to be physically based simplifications that allow us to solve linear PDE's or we will be a mathematical technique (e.g., a perturbation method) that linearizes the equations.  In the present problem the vanishing inertia, which allows us to neglect the left side of the equations results in a linear problem.  

2.  Solving PDE's:  You may not have solved a set of coupled partial differential equations before.  It is not as hard as it sounds.

	i.  The general procedure is to turn a PDE into an ordinary differential equation -- preferably an equation that you know how to solve.  In this case (and in many situations) this is done by Separation of Variables.  Separation of variables involves using a solution form, say F(x,y)  that is a function which is product of the functions of each of the independent variables, say X(x) Y(y).
	ii.  This assumed form of the solution is substituted into the equations and if it is going to work, you will either have one of functions appearing as an identical factor in every term (the case for this problem) or you will be able to group all of the X(x) functions together and all of the Y(y) terms together.  Since x and y can be varied independently, each of these groups of terms should be equation to the same constant.   
	
3.  You may not have solved a set of coupled partial ordinary equations before.  It is not as hard as it sounds.
	
	i.  You solve these by eliminating dependent variables (just as you would for a set of algebraic equations) and derivatives of independent variables in terms of the variable that you have chosen.
	ii.   If a dependent variable, that you need to eliminate (in this case Pressure) appears only as a derivative in two similar complex equations, it is often possible to take an extra derivative in one or both of the equations to create an identical term in each equation.  Then the equations can be subtracted to eliminate the variable.    

4.  Once you have gone through the process, you will have created an Euler differential equation.

5.  An Euler equation can be solved by assuming that the solution is a polynomial in the independent variable.

Mathematical Formulation

Clear[R,rr]

In case you ran the Peng-Robinson program without restarting.

Governing equations

Here are the equations from Truskey at al., p 123 and 134.  Note that I use “u” instead of “v” in the analysis below.

The momentum equation for velocity in the radial direction, (Subscript[u, r]) is:

ureq = D[p[r, θ], {r}] + μ*(-(-((2*Cot[θ]*uθ[r, θ])/r^2) - (2*D[uθ[r, θ], {θ}])/r^2 - 
      (2*ur[r, θ])/r^2 + D[r^2*D[ur[r, θ], {r}], {r}]/r^2 + 
      D[Sin[θ]*D[ur[r, θ], {θ}], {θ}]/(r^2*Sin[θ]))) + 
   ρ*(-(uθ[r, θ]^2/r) + ur[r, θ]*D[ur[r, θ], {r}] + (uθ[r, θ]*D[ur[r, θ], {θ}])/r)

The momentum equation for velocity in the tangent direction, (Subscript[u, θ]) equation is:

utheq = D[p[r, θ], {θ}]/r - (μ*(D[r^2*D[uθ[r, θ], {r}], {r}] + 
      2*D[ur[r, θ], {θ}] + D[D[Sin[θ]*uθ[r, θ], {θ}]/Sin[θ], 
       {θ}]))/r^2 + ρ*((uθ[r, θ]*ur[r, θ])/r + 
     ur[r, θ]*D[uθ[r, θ], {r}] + 
     (uθ[r, θ]*D[uθ[r, θ], {θ}])/r)

The mass balance (continuity) equation is

conteq=D[r^2*ur[r, θ], {r}]/r^2 + D[Sin[θ]*uθ[r, θ], {θ}]/
    (r*Sin[θ])

Dimensionless equations

How to make differential equations dimensionless

Somewhat more detail about how to make differential equations dimensionless using the chain rule is available in the notebook, Making a differential equation dimensionless.

This is done by using the chain rule for derivatives 


(∂/(∂r))= (ds/dr)(∂/(∂s))= (1/R)(∂/(∂s))
((∂^2)/(((∂r))^2))= ((((ds/dr)))^2)((∂^2)/(((∂s))^2))= (1/(R^2))((∂^2)/(((∂s))^2))


and making the substitutions:

s= r/R,
(Subscript[(OverTilde[u]), θ])= (Subscript[u, θ])/U, 
(Subscript[(OverTilde[u]), r])= (Subscript[u, r])/U

ureq

ureqx=ureq/.{
ur[r,θ]->(OverTilde[u])r[s,θ] U,uθ[r,θ]->(OverTilde[u])θ[s,θ] U,r->R s}

ureq1 = ureq /. {Derivative[a1_, a2_][ur][r, θ] :> (U*D[OverTilde[ur][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][uθ][r, θ] :> (U*D[OverTilde[uθ][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][p][r, θ] :> D[p[s, θ], {s, a1}, {θ, a2}]/R^a1, ur[r, θ] -> U*OverTilde[ur][s, θ], 
    uθ[r, θ] -> U*OverTilde[uθ][s, θ], r -> R*s}

First do the (Subscript[u, r]) equation

ureq1 = ureq /. {Derivative[a1_, a2_][ur][r, θ] :> (U*D[OverTilde[ur][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][uθ][r, θ] :> (U*D[OverTilde[uθ][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][p][r, θ] :> D[p[s, θ], {s, a1}, {θ, a2}]/R^a1, ur[r, θ] -> U*OverTilde[ur][s, θ], 
    uθ[r, θ] -> U*OverTilde[uθ][s, θ], r -> R*s}

ureq2=Apart[ureq1 R^2/μ /U]

ureq3=Collect[ureq2,ρ]

ureq4 = ureq3 /. Derivative[a1_, a2_][p][s, θ] :> D[OverTilde[p][s, θ], {s, a1}, {θ, a2}]*μ*(U/R)

It is seen that at this point, the factor ρ R U/μ, which is the Reynolds number,  is present in  three terms.  None of the other terms have any of these variables.  Since the equation is dimensionless, all of the terms have a magnitude of about one, or are much less than this or exactly zero.  Thus the value of the Reynolds number (in this case much less than one) determines which terms should be be kept for the solution.  

For the present problem we will eliminate the terms containing the Reynolds number by taking the limit of μ->∞.  This is the same as taking the Reynolds number going to 0.

back to objectives(physical)

back to conclusions

ureq5=Limit[ureq4,μ->Infinity]

ureq6=Simplify[ureq5]

Next do the (Subscript[u, θ]) equation

utheq1=utheq/.{Derivative[a1_, a2_][ur][r, θ] :> (U*D[OverTilde[ur][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][uθ][r, θ] :> (U*D[OverTilde[uθ][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][p][r, θ] :> D[p[s, θ], {s, a1}, {θ, a2}]/R^a1, ur[r, θ] -> U*OverTilde[ur][s, θ], 
    uθ[r, θ] -> U*OverTilde[uθ][s, θ], r -> R*s}

utheq2=Apart[utheq1 R^2/μ /U]

utheq3=Collect[utheq2,ρ]

utheq4=utheq3/.Derivative[a1_,a2_][p][s,θ]:>(∂_({s,a1},{θ,a2}))(OverTilde[p])[s,θ]μ U/R

It is seen that at this point, the factor ρ R U/μ, which is the Reynolds number,  is present in  three terms.  None of the other terms have any of these variables.  Since the equation is dimensionless, all of the terms have a magnitude of about one, or are much less than this or exactly zero.  Thus the value of the Reynolds number (in this case much less than one) determines which terms should be be kept for the solution.  

For the present problem we will eliminate the terms containing the Reynolds number by taking the limit of μ->∞.  This is the same as taking the Reynolds number going to 0.

utheq5=Limit[utheq4,μ->Infinity]

utheq6=Simplify[utheq5]

Make the continuity equation dimensionless

conteq1=conteq/.{Derivative[a1_, a2_][ur][r, θ] :> (U*D[OverTilde[ur][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][uθ][r, θ] :> (U*D[OverTilde[uθ][s, θ], {s, a1}, {θ, a2}])/R^a1, 
    Derivative[a1_, a2_][p][r, θ] :> D[p[s, θ], {s, a1}, {θ, a2}]/R^a1, ur[r, θ] -> U*OverTilde[ur][s, θ], 
    uθ[r, θ] -> U*OverTilde[uθ][s, θ], r -> R*s}

conteq2=Apart[conteq1 R /U]

The equations to be solved

ureq6

utheq6

conteq2

Solution of the equations with the boundary conditions

Solution of differential equations

Solution procedure for three coupled partial differential equations

We have three linear equations and three unknowns, the question is, how to solve this system of pdes.  The obvious answer is that we need to reduce them to ODE's that we know how to solve.  OK, how do we do this?  

We might expect that the general form of the solution is S(s) Θ(θ).  If this starts to work there are formal expansion techniques in terms of orthogonal polynomials that could allow us to proceed with no other knowledge of the solution.  However, this would be extremely tedious and involve use of mathematics that you probably don't know.  

Thus we invoke physics.  

In the low Reynolds number limit, the qualitative behavior of the flow field is very much determined by the boundaries because the fluid has no inertia and thus flows only where pushed.  In this case the "boundary" is that far away the flow is in a straight line.  Thus we look at the flow field at infinity to see its functional dependence.  It could be that this will provide insight into the shape.  

From the far field solutions, the likely form of the solution is

usr=Cos[θ] F[s];

usθ=Sin[θ] G[s];

back to objectives(physical)

back to objectives(mathematical)

back to conclusions

Go through the equations using the substitution to make ODE's

Start with the continuity equation.  We make the substitutions and we are pleased to see that the θ dependance of every term is the same.  So far the "guessed" form of the solution is working.

conteqsoln=conteq2/.{Derivative[a1_,a2_][OverTilde[ur]][s,θ]:>D[usr,{s,a1},{θ,a2}],Derivative[a1_,a2_][OverTilde[uθ]][s,θ]:>D[usθ,{s,a1},{θ,a2}],OverTilde[ur][s,θ]->usr,OverTilde[uθ][s,θ]->usθ}

back to objectives(mathematical)

We will need to substitute for G[s] in terms of F[s] in the momentum equations so we can end up with a single dependent variable, F[s] in the final equation.

geeofs=Simplify[G[s]/.(Solve[conteqsoln==0,G[s]])[[1]]]

Now let's work on the momentum equations,

ursoln1=ureq6/.{Derivative[a1_,a2_][OverTilde[ur]][s,θ]:>D[usr,{s,a1},{θ,a2}],Derivative[a1_,a2_][OverTilde[uθ]][s,θ]:>D[usθ,{s,a1},{θ,a2}],OverTilde[ur][s,θ]->usr,OverTilde[uθ][s,θ]->usθ}

ursoln2=ursoln1/.{G[s]->geeofs,Derivative[a1_][G][s]:>(∂_({s,a1}))geeofs}

uthsoln1=utheq6/.{Derivative[a1_,a2_][OverTilde[ur]][s,θ]:>D[usr,{s,a1},{θ,a2}],Derivative[a1_,a2_][OverTilde[uθ]][s,θ]:>D[usθ,{s,a1},{θ,a2}],OverTilde[ur][s,θ]->usr,OverTilde[uθ][s,θ]->usθ}

back to objectives(mathematical)

We can now eliminate G(s) and its derivatives with the results from the continuity equation.

uthsoln2=uthsoln1/.{G[s]->geeofs,Derivative[a1_][G][s]:>(∂_({s,a1}))geeofs}

Eliminate pressure leaving a single dependent variable, F[s].  This increases the degree of the derivatives as a trade.

Now we need to combine these equations to eliminate something.  The obvious thing is pressure.  The way to do it is to take the θ derivative of the r equation and the r derivative of the the θ equation.  Then subtract them to eliminate what now is the same 2nd derivative term of pressure.

back to objectives(mathematical)

ursoln3=D[ ExpandAll[ ursoln2],θ]

uthsoln3=D[FullSimplify[ s  uthsoln2],s]

Solve the single fourth order Euler equation.

When this is done, the result is an Euler type differential equation.  An Euler equation has the property that the power of the dependent variable multiplying a derivative, matches the power of the derivative.

Feq1=ExpandAll[-2 s^2 (ursoln3-uthsoln3)/Sin[θ]]

back to objectives(mathematical)

This solution can be obtained easily by assuming that F(s) = (s^n)

The solution, which matches Middleman, is:

Feq1

Fofs=F[s]/.DSolve[Feq1==0,F[s],s][[1]]

Boundary conditions in terms of F

Radial velocities

The 0 radial velocity on the surface, s=1, is F=0.

fsbceq1=Fofs/.s->1

To match the free stream radial velocity at s-->∞,F=1.  We can't directly substitute in s=∞, because of the s2 term.  Actually, (Subscript[c, 3]) will have to be 0.  Thus as s-->∞, F[s] becomes just (Subscript[c, 4]).

fsbceq2= C[4]

Tangential velocities

The 0 tangential velocity on the sphere is G=0 on s=1.  Matching the free stream far away is G=-1 for s->∞.  

In terms of F[s] these are:
F'[s]=0 @ s=1,
F'[s]=1 @ s->∞.

fsbceq3=D[Fofs,s]/.s->1

We don't do anything about the other infinity boundary condition because we already have (Subscript[c, 3])=0 there, so that the derivative of (Subscript[c, 3])is automatically 0.  The two far away boundary conditions are fit with the same constant because of the applicability of separation of variables solution and the angular functional forms.

Use the boundary equations to solve for the constants in the solution of the flow field.

We can solve the three equations to get the three constants.

cees=Solve[{fsbceq1==0, fsbceq2==1,fsbceq3==0},{C[1],C[2],C[4]}]/.C[3]->0

fsans=(Fofs/.cees[[1]])/.C[3]->0

gsans=ExpandAll[((geeofs/.{F[s]->fsans,F'[s]->D[fsans,s]})/.cees[[1]])/.C[3]->0]

Calculation of the pressure field

Both of these, of course, match Middleman.  We now need to get the pressure field.  It would have been easy to have obtained the angular form from the (Subscript[u, r]) equation just by substituting in the angular expressions for (Subscript[u, r]) and (Subscript[u, θ]).  However, we did not do this, so we go back to the last time we saw the pressure in the momentum equations.  This was the θ momentum equation.

Simplify[uthsoln1/Sin[θ]]

We can then substitute for G[s] and F[s],

psoln1=Simplify[uthsoln1/.{F[s]->fsans,F'[s]->D[fsans,s],G[s]->gsans, G'[s]->D[gsans,s],
G''[s]->D[gsans,{s,2}]}]

It should be apparent that the angular dependence is cos(θ) (because we take a θ derivative of p in this equation).  We need the solution to match a boundary condition on pressure (recall that pressure appears as first derivatives).   The obvious one for the r direction is that the pressure much match the far away pressure for the undisturbed flow, (Subscript[(OverTilde[p]), ∞]).

p2soln=(OverTilde[p])[s,θ]/.DSolve[ Expand[s psoln1]==0,(OverTilde[p])[s,θ],{s,θ}][[1]]

p3soln=p2soln/.C[1][s]->(Subscript[(OverTilde[p]), ∞])

We can transform this back to dimensional pressure for later calculation.

pfinalsoln=Expand[(p3soln/.{s->r/R}) μ U/R]

Velocity profiles

Let's calculate the velocity profiles from the answers that we have.

ur1=usr/.F[s]->fsans

urfinal=U ur1/.s->r/R

uth1=usθ/.G[s]->gsans

uthfinal=U uth1/.s->r/R

Examination of the solution

Contour plot of the velocity field

At this point it is interesting to look at the velocity field.  A way to do it is to look at the mean velocity as a function of the distance from the sphere and the angular location.  This is done by computing a total velocity, (vr^2) + ((vθ)^2).

ufield1 =Sqrt[(urfinal^2)+(uthfinal^2)]

Use the polar to Cartesian transformation:

ucart=ufield1/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[(y/x)]}

plot1=ContourPlot[ucart/. {R->1,U->1},{x,-5,5},{y,-5,5},Contours->12]

We can check the value of the velocity anywhere in the flow field.  We see that there is nowhere (outside of the sphere at least!), where the value exceeds 1.

Ω=Disk[{0,0},{1,1}];

sphere =Graphics[Ω]

N[ucart/.{x->.01,y->1}]/.{U->1,R->1}

N[ucart/.{x->.01,y->10}]/.{U->1,R->1}

N[ucart/.{x->-1,y->2}]/.{U->1,R->1}

The result is shown here with the sphere blacked out (and unfortunately some frame).  Note that the contours smoothly change as we approach the sphere.  The darker, the shading, the lower the velocity.  The point that could be surprising is that the fluid close to the sphere is always slower than far away.  There is no region where the fluid "speeds up" to get around the sphere which might happen for your "intuitively flowing" fluid.

back to objectives(physical)

back to conclusions

Show[plot1,sphere]

Drag on the sphere

We now get to the point of the problem, find the drag on the sphere.

Form drag

We need to integrate the normal stress, resolved onto the flow direction all over the sphere.  The normal stress is -p + 2 μ ((∂(Subscript[u, r])[r])/(∂r)).  It must be resolved onto the flow direction which takes Cos[θ].  The area element for a sphere is (R^2)Sin[θ] dθ dϕ.  The dϕ provides a 2 π and the integral becomes:

2 π Integrate[-p + 2 μ D[u[r],r]Cos[θ] (R^2) Sin[θ] ,{θ,0,π}].

argument=((-p+2 μ D[u[r],r])Cos[θ] (R^2) Sin[θ])/.{p->pfinalsoln, D[u[r],r]-> D[urfinal,r]}

FullSimplify[argument/.r->R]

The result of this integration is the "form drag", that is drag caused by normal stress.

back to objectives(physical)

back to conclusions

formdrag=2 π Integrate[argument,{θ,0,π}]/.r->R

Skin drag

Now we need to integrate the shear stress, resolved onto the flow direction all over the sphere.  The shear stress is (Subscript[T, rθ])= -μ( r (∂/(∂r))(((Subscript[u, θ])/r))+ (1/r)((∂(Subscript[u, r]))/(∂r))).  It must be resolved onto the flow direction which takes Sin[θ].  The area element for a sphere is (R^2)Sin[θ] dθ dϕ.  The dϕ provides a 2 π and the integral becomes:

2 π Integrate[-p + 2 μ D[u[r],r]Cos[θ] (R^2) Sin[θ] ,{θ,0,π}].

( -μ( r (∂_r)((((Subscript[u, θ])[r,θ])/r))+ (1/r)(∂_θ)(Subscript[u, r])[r,θ])Sin[θ] (R^2) Sin[θ])

argumenttang=( -μ( r (∂_r)((((Subscript[u, θ])[r,θ])/r))+ (1/r)(∂_θ)(Subscript[u, r])[r,θ])Sin[θ] (R^2) Sin[θ])/.{ D[(Subscript[u, r])[r,θ],θ]-> D[urfinal,θ],(∂_r) (Subscript[u, θ])[r,θ]->D[uthfinal,r],(Subscript[u, θ])[r,θ]->uthfinal}

FullSimplify[argumenttang/.r->R]

The result of this integration is the "skin drag", that is drag caused by tangential stress.

back to objectives(physical)

back to conclusions

skindrag=2 π Integrate[argumenttang,{θ,0,π}]/.r->R

totaldrag=formdrag+skindrag

Contour plot of the pressure field

At this point it is interesting to look at the pressure field.

pfield =pfinalsoln

-((3 R U μ Cos[θ])/(2 (r^2)))+((U μ (Subscript[(OverTilde[p]), ∞]))/R)

pfield =%

Plot[%24/.{R->1, μ->1, U->1,θ->3 π/4, (Subscript[(OverTilde[p]), ∞])->10},{r,10,1},PlotRange->All]

Plot[%24/.{R->1, μ->1, U->1,θ->π/4, (Subscript[(OverTilde[p]), ∞])->10},{r,1,10},PlotRange->All]

Use the polar to Cartesian transformation:

pcart=pfield/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[(( y)/x)]}

pcart/.y->0

pplot=ContourPlot[pcart/. {R->1,U->1,μ->0.01`,(Subscript[(OverTilde[p]), ∞])->10},{x,-5,5},{y,-5,5},Contours->12, ContourShading->False]

sphere=ContourPlot[(x^2)+(y^2),{x,-1,1},{y,-1,1},Contours->1, DisplayFunction->Identity];

The result is shown here with the sphere blacked out (and unfortunately some frame).  Note that the contours smoothly change as we approach the sphere.  However, the plot could be misleading.  On the upstream side, the pressure is increasing and the highest pressure is at the stagnation point (r =R, θ=π).  On the downstream side the pressure is lowest at r = R, θ = 0 and then increases back to the far away value.  

The symmetric color scheme does not convey this.

back to objectives(physical)

back to conclusions

Show[pplot,sphere]

Compare the solution for Creeping flow to a very high Reynolds number flow

To contrast the solution for no inertia, we write down the average velocity for the inviscid (meaning viscous terms are not considered) flow case from a standard reference

vinvis=(((Cos[θ] (1-(1/(r^3)))))^2)+(((-(1/2) Sin[θ] (2+(1/(r^3)))))^2);

vcartinvis=vinvis/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[(y/x)]};

plot2=ContourPlot[vcartinvis,{x,-5,5},{y,-5,5}, DisplayFunction->Identity, Contours->15];

Show[plot2,sphere,DisplayFunction->$DisplayFunction]

The code is as the shading gets darker, the velocity is slower.  (The plot is messed up just outside the sphere. We see that coming in the y=0 axis from x=-Infinity, the velocity is decreasing to agree with the: no flow through the surface boundary condition.  Looking at the x=0 axis, coming in from y=∞, the velocity is increasing reflecting the deflection of fluid around the sphere.  This results in a Bernoulli effect with local fluid acceleration.  Below we check some numbers to verify these statements.

N[vcartinvis/.{x->.01,y->1}]

N[vcartinvis/.{x->-2,y->-.001}]

N[vcartinvis/.{x->100,y->100}]

The same question can be asked about the pressure field.  We will find that the pressure is always at or below the far away value for Re=0.  Substitute for the pressure field equation,

Pressure field for high Reynolds number:

In contrast, if a high velocity fluid is impinging on a surface, there will be a pressure increase.  For an inviscid flow, from Bernoulli's equation, we find that dp ~ (U0^2-u^2)/2, where we have already plotted u^2

pinvis=((1-vcartinvis)/2)

plot4=ContourPlot[pinvis,{x,-5,5},{y,-5,5},Contours->25]

Show[plot4,sphere]

So that the pressure is highest at y=0, x=1,-1 (i.e. the stagnation point) and decreases out to infinity.  The pressure is lowest, consistent with the Bernoulli effect, at y=1,-1,x=0, and increases out to y=infinity.

A stream function, which is a function that is always tangent to the velocity field, is another way to view the solution of this problem.  The Re = 0 solution can be used to write the stream function which is, -1/2f[r] (r^2) ((sin[θ])^2).

re0stream=(1/2) fsans (r^2) ((Sin[θ])^2)/.{s->r/R}

We can make a plot of this if we transform from polar to Cartesian coordinates

stream1=re0stream/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[(y/x)],R->1}

plot5=ContourPlot[stream1,{x,-5,5},{y,-5,5},Contours->15,ContourShading->False,DisplayFunction->Identity];

Show[plot5,sphere,DisplayFunction->$DisplayFunction]

Note that at y=4, the contours are still deformed, the effect of the sphere is felt very far away!  

Now compare the inviscid flow case:

invisstream=(1/2) ((Sin[θ])^2) ((r^2)-(1/r))/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[(y/x)]};

plot6=ContourPlot[invisstream,{x,-5,5},{y,-5,5},ContourShading->False,ContourStyle->{Dashing[{.02,.02}]},DisplayFunction->Identity,Contours->12]

Show[plot5,plot6,sphere,DisplayFunction->$DisplayFunction]

back to objectives(physical)

back to conclusions

This last plot is a bit messy but shows that inviscid flow is disturbed much less a few diameters from the sphere than Stokes flow.

Conclusions

1.  Creeping flow is a term used for a flow that have effectively no inertia.  In this case the inertia terms are neglected  and the solution is obtained from the resulting linear equations.  The Reynolds number is very much smaller than unity.
        
2.  The solution technique involves using a solution form that is deduced from boundaries of the flow field, faraway from the sphere. 

3.  Because viscous forces dominate the flow field, the fluid can never accelerate above the free stream value even if an obstacle causes the fluid to be squeezed.  Thus the velocity in the region of the sphere just slows down and then returns to the free stream value.  

4.  Both normal stresses and tangential stresses contribute to the drag on the sphere.  These can be termed form drag and skin drag.  Note that both of these are linear in the velocity (consistent with the linear governing equations) and fluid viscosity.  The density, and thus the Reynolds number, does not appear.       

5.  Consistent with the fluid not accelerating, the pressure never increases above the free stream value.  The fluid has no inertia that would cause a pressure increase as the fluid slows down.

6.  The velocity decays slowly (as (1/r)) and thus the disturbance is felt very far away from the sphere.  This makes it difficult to do a real experiment, in a reasonable size container, that allows that sphere to fall at a speed specified by the drag that is predicted from the analysis here.  The very high Reynolds number case decays much faster.