Creeping flow caused by a rotating 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

Mark.J.McCready.1@nd.edu
http://www.nd.edu/~mjm/


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, if you use it, this notice remain visible to other users. 
There is no charge for copying and dissemination.   

Version:  8/23/18
More recent versions of this notebook may be available at the web site:
http://www.chemeprof.com/cbe60544.html

Summary

This notebook shows how to solve the flow field for a rotating sphere at low Reynolds number showing the 0 order and first order solutions.  A regular perturbation expansion is used.  At the lowest order, inertial terms are neglected and the flow is pure rotation that decays away from the sphere.  The perturbation improvement takes into account weak inertia and this causes a secondary flow in the r and θ directions. 

This notebook is intended as a supplement to 
L. G. Leal  (1992) Laminar flow and Convective Transport Processes, Butterworth.
and
W. M. Deen (2012)  Analysis of Transport Phenomena, Oxford.

The references for this problem in the literature are:

Keywords

Creeping flow, Rotating Sphere, Regular Perturbation, Low Reynolds number fluid dynamics

Problem formulation

Consider an infinite quiescent fluid.  A sphere of radius a is placed in the region and caused to rotate with angular velocity, Ω, such that the Reynolds number, R = (a^2)Ω ρ/μ is small compared to 1.  The 0 order governing equation is just the (Subscript[v, ϕ]) equation without any inertia terms.  The answer will have the  same angular dependence as the velocity on the surface of the sphere (See Leal's discussion starting on Page 135, Deen: 321) and decay with distance.  Since inertia is weak, any radial flow associated with centrifugal effects will enter at higher order.

We solve the 0 order problem first.  Then we look to find out how the 0 order solution produces weak inertia to cause a flow at first order.

The 0 order problem.  Flow is only in ϕ direction.

A sphere is rotating in an otherwise quiescent medium.  The equation for (Subscript[v, ϕ]), assuming creeping flow, is:

vϕeq=(((((∂_r)((r^2) (∂_r)vϕ[r,θ]))/(r^2))+(((∂_θ)(((∂_θ)(vϕ[r,θ] Sin[θ]))/(Sin[θ])))/(r^2)))/R)

How to solve a PDE.

1.  Turn it into an ODE

If we are to obtain an analytic solution to a PDE, we will need to convert it to an ODE.  Two general methods are available.  

1.  Separation of variables, where the assumed form of a solution for, say, v(r,θ) is R(r) Θ(θ).  That is the solution is the product of two functions that are only a function of 1 independent variable.

This solution is then substituted into the PDE and rearranging is done until you get all of the R(r)'s on one side and all of the Θ(θ)'s on the other side.  If this can be done, each side is equal to a constant (the same one) and you can solve the two ODE problems independently.  Also, (most of the time) the equation that is being solved on a finite domain, will be  Sturm-Liouville eigen value problem.  After each ODE is solved with the boundary conditions, the final solution is obtained from the product of the two parts.

If you cannot get all of the r's and all of the θ's separated, then you are stuck.  Either the problem is not separable, you have made a mistake or the form of the trial solution that have assumed is not correct.  

2.  Combination of variables or use of a similarity variable.  If the problem has no finite domain for either independent variable and there is consequently no geometric length scale to use to nondimensionalize the problem, you expect that the solutions will exhibit self-similarity.  If this is the case, the PDE can often be converted to an ODE by a change of independent variable that is a dimensionless function of the two independent variables.  You can see this done in excruciating detail in http://www.nd.edu/~mjm/thermal_boundarylayer.nb.  

You don't need to worry about it right now however.  We can solve the present problem with Separation of Variables.

2.  How to get a functional form for one of the parts of the separation of variables solution to make the problem simpler.

Please take a quick peak here.

It turns out that with Mathematica's help, you can use a general separation of variables approach to this problem.

Here is the equation

vϕeq=(((((∂_r)((r^2) (∂_r)vϕ[r,θ]))/(r^2))+(((∂_θ)(((∂_θ)(vϕ[r,θ] Sin[θ]))/(Sin[θ])))/(r^2)))/R)

We assume the appropriate functional form

vfix=f[r] Q[θ];

vϕeqx1=vϕeq/.{Derivative[a1_,a2_][vϕ][r,θ]:>(∂_({r,a1},{θ,a2}))vfi,vϕ[r,θ]->vfi}

vϕeqx2= FullSimplify[vϕeqx1]

vϕeqx3 = ExpandAll[vϕeqx2/Q[θ]/f[r] r^2 R]

So that the first two terms equal a constant as do the last three terms.  Thus the equation is separable.  However, have fun trying to solve it.

A general principle for Creeping flows is that the flow field is largely determined by the boundaries. (See Leal).  Thus we expect to obtain a form of the separation of variables solution from the boundary conditions.  The key issue here is the θ dependence of velocity.  We will assume a sine dependence which is the same as velocity on the surface of the sphere.

return to conclusions

vfi=f[r] Sin[θ];

Now substitute and see what happens

vϕeqO01=vϕeq/.{Derivative[a1_,a2_][vϕ][r,θ]:>(∂_({r,a1},{θ,a2}))vfi,vϕ[r,θ]->vfi}

The substitution has worked and the angular dependence of the equation has simplified.

vϕeqO02=Expand[vϕeqO01]

We can put this into a standard form, an Euler differential equation, by multiplying by the appropriate factor

vϕeqO03=Expand[((vϕeqO02 (r^2) R)/(Sin[θ]))]

This is readily solved.

ansO01=DSolve[{vϕeqO03==0,f[a] == a Ω},f[r],r]

ansO01=DSolve[{vϕeqO03==0},f[r],r]

{{f[r]->r C[1]+(( (a^3)Ω)/(r^2))}}

Because we need the solution to decay away from the sphere and be 0 sufficiently far away, C[1] = 0.  On the surface of the sphere, the velocity is  a Ω Sin[θ].  The Sin[θ] is already built in so we just have that C[2] = (a^3) Ω

ansO02=f[r]/.ansO01[[1]]

ansO03=ansO02/.{C[1]->0}

(( (a^3)Ω)/(r^2))

This gives a solution where the velocity is only in the ϕ direction and decay as 1/(r^2).   The θ dependence is the same as the velocity on the sphere.

vphi0 =(((a^3) Ω)/(r^2)) Sin[θ]

Find the torque.  We need the stress component (Subscript[τ, rϕ]) which is:

stress = μ r D[vphi0/r,r]

Integrate[Integrate [stress  a Sin[θ] a^2 Sin[θ] ,{θ,0,π}],{ϕ,0, 2 π}]/.r->a

return to conclusions

The 1/(r^2) dependence is fast enough that a regular perturbation method can be used to capture the first effects of inertia.  

We should address the claim made by Deen on page 321 that it is not possible to find a pressure for the pure rotation solution because the two pressure terms cannot be reconciled.  (eqs 8.3-12 and 8.3-13).  I don’t think he is correctly ordering the various terms in the r and θ equations.  We are solving for the flow in the limit of 0 Re.  This is the only justification for being able to separately solve the vϕ equation by itself.  Following this ordering, if we examine the r and θ equations, we should then solve ∂(Subscript[p, 0])/∂r =0, 1/r ∂(Subscript[p, 0])/∂θ=0.  It is not correct to write his eq. 8.3-10 with the centrifugal term being of the same order.  

Thus we can proceed with a perturbation expansion and get a valid correction to the pure rotational flow!

Ratio of inertial and viscous terms

The viscous term is:

μ(((∂_r)((r^2) (∂_r)vϕ[r,θ]))/(r^2))

viscousterm = %/.{Derivative[a1_,a2_][vϕ][r,θ]:>(∂_({r,a1},{θ,a2}))vphi0,vϕ[r,θ]->vphi0}

An inertia term is

ρ (((vϕ[r,θ])^2)/r)

inertiaterm= %/.{Derivative[a1_,a2_][vϕ][r,θ]:>(∂_({r,a1},{θ,a2}))vphi0,vϕ[r,θ]->vphi0}

The ratio is

(inertiaterm/viscousterm)/.Ω-> μ Reω/a^2/ρ

So we see that as r/a-> Infinity, this ratio --> 0   (Uniformly valid!)

The first order problem which comes from perturbation theory

Now suppose that we want to improve this by using a perturbation solution.  If we examine the inertia (nonlinear) terms in the (Subscript[v, ϕ]) equation, we see that the 0 order solution for (Subscript[v, ϕ]) produces no correction terms at (R^1) order.  The pure rotation solution that was obtained is an exact solution to the (Subscript[v, ϕ])equation in the absence of (Subscript[v, r ])or (Subscript[v, θ]).  Thus we need to examine the (Subscript[v, r]) and (Subscript[v, θ]) equations to find the first improvement term.  In the (Subscript[v, r]) equation, there is a  -(Subscript[v, ϕ])2/r term.  In the θ equation, there is a similar term.  We are seeing that rotation with finite R causes a weak centrifugal force that leads to a flow in the r direction -- we suspect outward at the equator.  To conserve mass, there is necessarily an inflow someplace, presumably at the point of weakest inertia, the poles.  This flow at (R^1)requires both the (Subscript[v, r]) and (Subscript[v, θ]) equations to describe it.  Let's see what we can do. 


Here is the axisymmetric (about the ϕ axis) (Subscript[v, r]) equation with the additional inertia term.  This is produced by the expansion of the equations at (R^1) order.

return to conclusions

vr1eq=-(((vϕ[r,θ])^2)/r)-((∂_r)(((∂_r)(vr[r,θ] (r^2)))/(r^2))+(((∂_θ)(Sin[θ] (∂_θ)vr[r,θ]))/((r^2) Sin[θ]))-((2 (∂_θ)(vθ[r,θ] Sin[θ]))/((r^2) Sin[θ]))-(∂_r)p[r,θ])

Here is the corresponding (Subscript[v, θ]) equation.

vθ1eq=-((((vϕ[r,θ])^2) Cot[θ])/r)-((((∂_r)((r^2) (∂_r)vθ[r,θ]))/(r^2))+(((∂_θ)(((∂_θ)(vθ[r,θ] Sin[θ]))/(Sin[θ])))/(r^2))+((2 (∂_θ)vr[r,θ])/(r^2))-(((∂_θ)p[r,θ])/r))

Now we can derive the equation for the stream function easily enough, or could we do something else? (See Leal for a general method).

Let's get the governing equation in terms of the stream function.  

First, eliminate the pressure with a cross derivative.

eqO11=Expand[(∂_θ)vr1eq-(∂_r)(r vθ1eq)]

From the continuity equation we can replace (Subscript[v, r]) and (Subscript[v, θ])using the definition of the stream function.

eqO12= eqO11/.{Derivative[a1_,a2_][vr][r,θ]:>(∂_({r,a1},{θ,a2}))(((∂_θ)ψ[r,θ])/((r^2) Sin[θ])),vr[r,θ]->(((∂_θ)ψ[r,θ])/((r^2) Sin[θ])),Derivative[a1_,a2_][vθ][r,θ]:>(∂_({r,a1},{θ,a2}))(-(((∂_r)ψ[r,θ])/(r Sin[θ]))),vθ[r,θ]->-(((∂_r)ψ[r,θ])/(r Sin[θ]))};

Short[eqO12]

We can make it look a lot better.

eqO13=Expand[eqO12 Sin[θ]]

Now substitute the 0 order solution for (Subscript[v, ϕ])

eqO14=eqO13/.{vϕ[r,θ]->vphi0,(∂_θ)vϕ[r,θ]:>(∂_θ)vphi0,(∂_r)vϕ[r,θ]:>(∂_r)vphi0}

Now we have a PDE that we need to solve.    The problem is, how?
One way is to guess the function form for the θ dependence.  The boundary relation for the 0 order problem is clear but does readily not give insight into the first order problem.  Of course we could just guess and try.  But what is this equation really?

It should be an axisymmetric creeping flow problem with an inhomogeneous term

(E^4) ψ = 6 (a^6) (Ω^2) Cos[θ] ((Sin[θ])^2)/(r^5).

We know the general solutions to this equation.  We would need a particular solution.  The boundary conditions for the rotating problem are all fit at 0 order so the homogenous solution and the particular solution must cancel on the boundary for this perturbation correction to work.  

Now that we have some familiarity, we could translate using η≡cos(θ) (see page 160 of Leal, 458 in new version) and look at the general solutions.  A second way would be to look at the η solutions in terms of Gegenbauer polynomials for the expected form of the θ dependence(see page 163 (461)).  What are the choices:
Cos[θ] + 1, which is not in the general solution.

The next one is

Simplify[(Cos[θ]^2-1)/2]

It does not look promising.  Continue on with the third,

Simplify[Cos[θ]^3-Cos[θ]]

This looks great!! It has the same form as the inhomogeneous term.  Try this form for the θ dependence.

eqO15=Expand[eqO14/.{Derivative[a1_,a2_][ψ][r,θ]:>(∂_({r,a1},{θ,a2}))(((Sin[θ])^2) Cos[θ] f1[r])}]

eqO16= Simplify[eqO15]

Now we easily eliminate the θ dependence.

eqO17= Expand[(((r^4) eqO16)/(Cos[θ] ((Sin[θ])^2)))]

The equation is Euler's equation which is easily solved.

ans011=DSolve[eqO17==0,f1[r],r]

ans012=f1[r]/.ans011[[1]]

As mentioned earlier, we need the (R^1) solution to be 0 on the boundaries.  As r -->∞ this requires that C[3], C[2] =0.

ans013=ans012/.{C[2]->0,C[3]->0}

Now fit the boundary conditions on the sphere

bc1=(ans013/r)/.r->a

bc2=-(((∂_r)ans013)/(r^2))/.r->a

Everything is 0 on the boundaries

constsO1=Solve[{bc1==0,bc2==0},{C[1],C[4]}]

Thus, the stream function for (R^1) order is

ansO1last=(ans013/.constsO1[[1]]) ((Sin[θ])^2) Cos[θ]

ansO1last=(-(1/8) (a^5) (Ω^2)-(((a^7) (Ω^2))/(8 (r^2)))+(((a^6) (Ω^2))/(4 r))) Cos[θ] ((Sin[θ])^2)

This matches Collins (1955).

Examine the first order solution graphically.

Make the requisite polar to Cartesian transformation,

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

Substitute in some numbers,

anstoplot=ansO1Cartesian/.{Ω->5,a->1}

The stream lines look like:     (Note that it could have been that we would need to plot only part of the domain given the polar to Cartesian transformation.

ContourPlot[anstoplot,{x,-10,10},{y,-10,10},Contours->20,PlotPoints->50]

return to conclusions

We see the circulation pattern.  Flows are strongest outward along the equator and this fluid is replaced by inflows along the poles.

The velocity components are

vr1zz =Simplify[ (((∂_θ)ansO1last)/((r^2) Sin[θ]))]

vr1 =FullSimplify[ vr1zz]

vθ1 = -(((∂_r)ansO1last)/(r Sin[θ]))

Now what is the pressure?

vr1eq=-(((vϕ[r,θ])^2)/r)-((∂_r)(((∂_r)(vr[r,θ] (r^2)))/(r^2))+(((∂_θ)(Sin[θ] (∂_θ)vr[r,θ]))/((r^2) Sin[θ]))-((2 (∂_θ)(vθ[r,θ] Sin[θ]))/((r^2) Sin[θ]))-(∂_r)p[r,θ])

vr1eq=-(((vϕ[r,θ])^2)/r)-((∂_r)(((∂_r)(vr[r,θ] (r^2)))/(r^2))+(((∂_θ)(Sin[θ] (∂_θ)vr[r,θ]))/((r^2) Sin[θ]))-((2 (∂_θ)(vθ[r,θ] Sin[θ]))/((r^2) Sin[θ]))-(∂_r)p[r,θ])

vr1eq/.{ vϕ[r,θ]->(((a^3) Ω)/(r^2)) Sin[θ], Derivative[a1_,a2_][vr][r,θ]:>  ((∂_({r,a1},{θ,a2}))(  vr1)) , Derivative[a1_,a2_][vθ][r,θ]:>  ((∂_({r,a1},{θ,a2}))(  vθ1)),vr[r,θ]->vr1,vθ[r,θ]->vθ1 }

FullSimplify[%]

DSolve[{D[-p[r],r] == (((a^5) (Ω^2) (-3 r+(16 a-9 r) Cos[2 θ]))/(8 (r^5)))},p[r],r]

pr1zz = (p[r]/.%[[1]])/.C[1]->0

Again a match to Collins.

pr1zz =-(((a^5) (Ω^2) (r+(-4 a+3 r) Cos[2 θ]))/(8 (r^4)))

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

Substitute in some numbers,

p01plot=p01Cartesian/.{ω->5,Ω->5,R->1, Re->.01,a->1}

The stream lines look like:     (Note that it could have been that we would need to plot only part of the domain given the polar to Cartesian transformation.

ContourPlot[p01plot,{x,-10,10},{y,-10,10},Contours->20,PlotPoints->50]

What is the pressure from the vθ equation?

vθ1eq=-((((vϕ[r,θ])^2) Cot[θ])/r)-((((∂_r)((r^2) (∂_r)vθ[r,θ]))/(r^2))+(((∂_θ)(((∂_θ)(vθ[r,θ] Sin[θ]))/(Sin[θ])))/(r^2))+((2 (∂_θ)vr[r,θ])/(r^2))-(((∂_θ)p[r,θ])/r))

vθ1eq=-((((vϕ[r,θ])^2) Cot[θ])/r)-((((∂_r)((r^2) (∂_r)vθ[r,θ]))/(r^2))+(((∂_θ)(((∂_θ)(vθ[r,θ] Sin[θ]))/(Sin[θ])))/(r^2))+((2 (∂_θ)vr[r,θ])/(r^2))-(((∂_θ)p[r,θ])/r))

vθ1eq/.{ vϕ[r,θ]->(((a^3) Ω)/(r^2)) Sin[θ], Derivative[a1_,a2_][vr][r,θ]:>  ((∂_({r,a1},{θ,a2}))(  vr1)) , Derivative[a1_,a2_][vθ][r,θ]:>  ((∂_({r,a1},{θ,a2}))(  vθ1)),vr[r,θ]->vr1,vθ[r,θ]->vθ1 }

FullSimplify[r %]

DSolve[{(((a^5) (4 a-3 r) (Ω^2) Sin[2 θ])/(4 (r^4)))+Derivative[0,1][p][r,θ]==0},p[r,θ],θ]

The constant is a function of r.  We need the result from the “r” equation to get it.  Let’s call it gg[r]

pr1θ = (((a^6) (Ω^2) Cos[2 θ])/(2 (r^4)))-((3 (a^5) (Ω^2) Cos[2 θ])/(8 (r^3))) + gg[r]

pr1θθ= (((a^6) (Ω^2) Cos[2 θ])/(2 (r^4)))-((3 (a^5) (Ω^2) Cos[2 θ])/(8 (r^3)))+gg[r]

Here is the pressure from the r equation:

pr1zz

Expand[pr1zz-pr1θθ]

So the pressure from the r and θ will match if we add the correct function of r as the constant from the θ integration, gg[r] =  -(((a^5) (Ω^2))/(8 (r^3)))

pr1θ=(((a^5) (4 a-3 r) (Ω^2) Cos[2 θ])/(8 (r^4)))-(((a^5) (Ω^2))/(8 (r^3)))

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

Substitute in some numbers,

p01plot=p01Cartesian/.{ω->5,Ω->5,R->1, Re->.01,a->1}

The stream lines look like:     (Note that it could have been that we would need to plot only part of the domain given the polar to Cartesian transformation.

ContourPlot[p01plot,{x,-10,10},{y,-10,10},Contours->20,PlotPoints->50]

Can we get a correction to vphi?

The (Subscript[v, ϕ]) pde at this order will have the inertia terms: (Subscript[v, ϕ0])*(Subscript[v, r1]), (Subscript[v, ϕ0])*(Subscript[v, θ1]) and (Subscript[v, ϕ2]) terms.

vϕeq2= vr[r,θ] D[vϕ[r,θ],r] + vθ[r,θ] /r D[vϕ[r,θ],θ] + vr[r,θ] vϕ[r,θ]/r + vθ[r,θ] vϕ[r,θ] Cot[θ]/r-(((∂_r)((r^2) (∂_r)vϕ2[r,θ]))/(r^2))-(((∂_θ)(((∂_θ)(vϕ2[r,θ] Sin[θ]))/(Sin[θ])))/(r^2))

vϕeq21=vϕeq2/.{vϕ[r,θ]->vphi0,(∂_θ)vϕ[r,θ]:>(∂_θ)vphi0,(∂_r)vϕ[r,θ]:>(∂_r)vphi0}

vϕeq22= vϕeq21/.{Derivative[a1_,a2_][vr][r,θ]:>(∂_({r,a1},{θ,a2}))(((∂_θ)ansO1last)/((r^2) Sin[θ])),vr[r,θ]->(((∂_θ)ansO1last)/((r^2) Sin[θ])),Derivative[a1_,a2_][vθ][r,θ]:>(∂_({r,a1},{θ,a2}))(-(((∂_r)ansO1last)/(r Sin[θ]))),vθ[r,θ]->-(((∂_r)ansO1last)/(r Sin[θ]))}

ExpandAll[%]

For vϕ, we expect the same basic angle function, Sin[θ].

vϕeqO01=vϕeq22/.{Derivative[a1_,a2_][vϕ2][r,θ]:>(∂_({r,a1},{θ,a2}))(Sin[θ]vf2[r]),vϕ2[r,θ]->Sin[θ]vf2[r]}

vϕeqO02 = Expand[vϕeqO01/Sin[θ]]

Since the 0 order solution fits the boundary condition perfectly, we need this order to be 0 on the sphere

vϕeqO03 = DSolve[{vϕeqO02==0,vf2[a] ==0,vf2[∞]==0},vf2[r],r]

vf2ans = vf2[r]/.%[[1]]

Don’t forget to put the angle function back in!

vphi2 = Sin[θ]ExpandAll[%]

The r-ϕ stress component to get the torque is:

stress2 = μ r D[vphi2 /r,r]

Do the appropriate integrations on the sphere...

Integrate[Integrate [stress2  a Sin[θ] a^2 Sin[θ] ,{θ,0,π}],{ϕ,0, 2 π}]/.r->a

The ‘150 π” looks promising....

N[%/.{a->1, Ω->1,μ->1}]

This matches the result from W. D. Collins, Mathematika 2 (1955) p 42-47!

π/150.

This matches the boundary conditions so maybe it is correct.

Let’s see what it looks like.

vfi2 = (vf2ans Sin[θ])

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

Substitute in some numbers,

anstoplot=ansO1Cartesian/.{Ω->5,a->1}

The stream lines look like,

ContourPlot[anstoplot,{x,0,10},{y,0,10},Contours->20,PlotPoints->50]

It is hard to tell if this is correct!  The lighter colors should be the O (R^2) enhancement to the 1/r^2 flow.

Plot[{vfi2/.{Ω->5,a->1, θ->π/2},vfi2/.{Ω->5,a->1, θ->π/3},vfi2/.{Ω->5,a->1, θ->π/4},vfi2/.{Ω->5,a->1, θ->π/8},vfi2/.{Ω->5,a->1, θ->0}},{r,1,10},PlotLegends->"Expressions"]

Let’s get the mass transfer enhancement

We realize that the solution should be of the form:

C[r,θ] = (Subscript[C, 0])[r] + Pe (Subscript[C, 1])[r,θ]

Here is the 0 order solution for pure diffusion

ϕ0diffusioneq= diff 1/r^2 D[r^2 D[c[r],r],r]

DSolve[{%==0},c[r],r]

(c[r]/.%[[1]])/.C[2]->0

Solve[(%/.r->1) ==c0,C[1]]

cpurediffusion = %%/.%[[1]]

Let’s define a Nusselt number using the standard convention.  Sh or Nu =2 for pure diffusion.

(-2 diff D[%316,r](4 π a^2)/(4 π a^2)/diff/(c0) )/.r->1

The order 1 solution with enhancement.  We “sort of” get an answer.

If the sphere is rotating slowly such that no secondary currents are created, mass transfer will be by pure diffusion.  This is because in the convection-diffusion equation vϕ shows up only in a term that is otherwise symmetric in ϕ.

Thus vr and vθ will contribute to mass transfer, but not vϕ.

Here is the order Re^1 solution, from above for vr and vθ.

streamfunctionRe1=(-(((a^7) (Ω^2))/(5 (r^2)))+(((a^6) (Ω^2))/(4 r))-(1/20) (a^2) (r^3) (Ω^2)) Cos[θ] ((Sin[θ])^2)

vr1=FullSimplify[ (((∂_θ)streamfunctionRe1)/((r^2) Sin[θ]))]

vθ1 =FullSimplify[ -(((∂_r)streamfunctionRe1)/(r Sin[θ]))]

To get the enhancement from the secondary flow, we substitute the R^1 solution into the left side for vr and vθ

Here is the pde

ϕdiffusioneq= vr[r,θ] D[C0[r,θ],r] + vθ[r,θ] /r D[C0[r,θ],θ]- diff(1/r^2 D[r^2 D[C1[r,θ],r],r]+
1/r^2/Sin[θ]D[Sin[θ] D[C1[r,θ],θ],θ])

Here we substitute the velocity components.

difeq1 = ϕdiffusioneq/.{vr[r,θ]->vr1, vθ[r,θ]->vθ1}

Do this to grab just the convection terms

temp1=Expand[difeq1]/.diff->0

Now we replace the concentration derivatives with the order 0 (i.e. Pe^0) solution for the concentration field.  This is just pure diffusion.

temp2 = temp1/.{Derivative[a1_,a2_][C0][r,θ]:>(∂_({r,a1},{θ,a2}))(c0/r),C0[r,θ]->(c0/r)}

Now restore the diffusion terms.

temp3 = temp2 -diff (((Csc[θ] (Cos[θ] Derivative[0,1][C1][r,θ]+Sin[θ] Derivative[0,2][C1][r,θ]))/(r^2))+((2 r Derivative[1,0][C1][r,θ]+(r^2) Derivative[2,0][C1][r,θ])/(r^2)))

We would like to solve this equation by turning it into an ODE.  We need an angle function.

First try is Sin[θ]

difeq2=temp3/.{Derivative[a1_,a2_][C1][r,θ]:>(∂_({r,a1},{θ,a2}))(Sin[θ] f2[r]),C1[r,θ]->(Sin[θ]f2[r])}

difeq3=ExpandAll[difeq2/Sin[θ]]

It does not appear that this operator will completely separate.  (If you find a function, please let me know!!).  Note the term: -((diff ((Cot[θ])^2) f2[r])/(r^2))

However, DSolve is pretty good so just let it go!

f2answer1 =DSolve[{difeq3==0, f2[a] ==0}, f2[r],r]

Wow!  I’ll bet there is a lot of secrets in this!

Since I know that the solution is 0 on the sphere, I am going to take C[1]=0 with the assumption that it will be multiplying the highest order term.  

We don’t need the diffusivity as it would have been scaled away in the perturbation expansion.

f2answer = (f2[r]/.f2answer1[[1]])/.{a->1,C[1]->0,diff->1}

Let’s get the flux

flux2 = D[f2answer Sin[θ],r]/.r->1

We could try to get an average, but after several minutes, I gave up.

Timing[fluxaverage= Integrate[flux2 Sin[θ],{θ,0,π/2}]]

ans = %23[[2]]

Let’s just do the series expansion....

seriesflux2 = Series[flux2 Sin[θ],{Ω,0,2}]

We punt and pick just a single angle.  We see that the first term is Ω^2.  This means that the correction is O Pe^2.

seriesflux2/.θ->π/2

And get a numerical value for the first term

N[%]

When I had fixed θ=π/2 before solving the pde I got”

((2 c0 (Ω^2))/(-15+(Sqrt[5])))+((O[Ω])^3)

N[%]

Which is only a little different.

I could try one more thing.  But this did not give an answer after several minutes.

Integrate[ Normal[seriesflux2] Sin[θ],{θ,0,π/2}]

Thus the Nusselt number goes as  2+ Pe^2.

For flow past a sphere the Nu = 2 + Pe 

So the streaming flow causes more enhancement.

Conclusions

1.  As expected for creeping flows, the form of the solution is determined by the boundary conditions (because the fluid has no inertia to create its own dynamics).  Thus the θ dependence of the 0 order solution is given by the change in the velocity with angle on the surface of the sphere.

2.  The zero order solution is that the fluid is dragged around by the rotating sphere and decays as 1/(r^2). 

3.  We might try to find a correction to the (Subscript[v, ϕ]) equation by expanding (Subscript[v, ϕ])in powers of R and substituting in the (Subscript[v, ϕ]) equation.  However, this does not produce a new problem to solve.  The reason is that the 0 order solution is an exact solution to the (Subscript[v, ϕ])equation. 

4.  The first effect of inertia is to cause a circulation in the r and θ directions.  You might expect that inertia is largest at the equator and thus this is where the outflow occurs.  To match the outflow, there is an inflow near the poles.