Creeping flow caused by a rotating sphere

This notebook has been written in Mathematica by 

Mark J. McCready
Professor and Chair 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:  6/8/99
More recent versions of this notebook should be available at the web site:
http://www.nd.edu/~mjm/rotatingsphere.nb

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.

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) 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)

We can assume a form of the solution from the boundary conditions.  Note that this is simply saying the velocity is highest at the equator and 0 at the "poles".

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 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[r], r]

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, 
          C[2] -> a^3  Ω}

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.

O0ans = ansO03 Sin[θ]

return to conclusions

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

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 stream function

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

From the continuity equation we can replace (Subscript[v, r]) and (Subscript[v, θ]).

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, θ] -> 
            O0ans, (∂_θ) vϕ[
                r, θ] :> (∂_θ) 
              O0ans, (∂_r) vϕ[
                r, θ] :> (∂_r) O0ans}

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 equation.  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) 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).  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]

By knowing what we want we can simplify the equation,

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[4] =0.

ans013 = ans012 /. {C[3] -> 0, C[4] -> 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[2]}]

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

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

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,

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.

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.