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:  4/1/17
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

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  R Sin[θ] R^2 Sin[θ] ,{θ,0,π}],{ϕ,0, 2 π}]/.r->R