Inviscid flow past a stationary 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, that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  6/18/00
More recent versions of this notebook should be available at the web site:
http://www.nd.edu/~mjm/InviscidSphere.nb

Reference:  M. M. Denn, Process Fluid Dynamics, Prentice Hall, 1978

Summary

This notebook considers the problem of an inviscid flow past a stationary sphere.  Solutions are obtained using both a potential function and a stream function.  The pressure field around the sphere and the drag, which is 0 (d'Alembert's paradox) are found.

This solution is a nice example of how for an inviscid flow (the idealized limit of a very high Reynolds number flow), the flow field can be computed even though the governing Euler equations are nonlinear.  This is because the flow can be considered irrotational (vorticity=0) so that the flow field can be defined as the gradient of a scalar potential function.  

Such a formulation should give useful results for high reynolds number flows away from boundaries where vorticity is produced.  Perhaps the best example of the utility of inviscid flow theory is for flow around airfoils to give the pressure field and thus the lift.  The drag is obtained from boundary-layer theory -- where viscous effects are important.

Mathematica aside

I have used a mix of input notation to show how it can be done.   Also, some of the major equations are shown in input form and also in typeset form for clarity.  The typeset input equations are usually marked as inactive and so are not read in by Mathematica .. but they could be.

I would enter the dynamic boundary condition from a key pad in input form as

\(dc = \(-\ γ\)\ D[η[x, t], {x, 
              2}]\  - 
		\(ρ\_1\) \((\(-D[ϕ\_1[x, y, t], t]\) - 
              U1\ D[ϕ\_1[x, y, t], x] - 
              g\ η[x, 
                  t])\) + 
			\(ρ\_2\) \((\(-D[ϕ\_2[x, y, t], 
                  t]\) - 
					U2\ D[ϕ\_2[x, y, t], x] - 
              g\ η[x, t])\)\)

In Standard form, which is a little shorter, you would need the type set window to make this practical.

\(dc = \(-\((γ\ ∂\_{x, 2}η[x, 
                  t])\)\) - ρ\_1\ \((\(-∂\_t ϕ\_1[x, y, 
                    t]\) - U1\ ∂\_x ϕ\_1[x, y, t] - 
              g\ η[x, 
                  t])\) + ρ\_2\ \((\(-∂\_t ϕ\_2[x, y, 
                    t]\) - U2\ ∂\_x ϕ\_2[x, y, t] - 
              g\ η[x, t])\)\)

Finally it may be more easily read in Traditional form as

dc=-(γ (((∂^2) η (x, 
                        t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x, y, 
                        t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x, y, 
                        t])/(∂x))-g η (x, t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x, y, 
                        t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x, y, 
                        t])/(∂x))-g η (x, t))

Note that you can convert any "cryptic" input expression by selecting the cell, going up to the cell menu and selecting "convert to traditional form".  This is "shift curly t" on the keypad.

Keywords

Inviscid flow, Ideal flow theory, Irrotational, Transport phenomena,  Fluid Dynamics

Homework questions for this notebook:

Run the file, read the exposition and then answer these questions.
 
1.	  For an inviscid flow, what physical observation is made that leads to the governing equation for the flow field in terms of a scalar function?

2.	Why is there a substantial simplification of the original problem in terms of the N-S equation with η =0., Re--> ∞?

3.	Flow separation occurs readily in the flow past a sphere problem if Re is sufficiently high.  Explain why it does not occur for an inviscid flow.

Problem statement

The problem of interest is a uniform, inviscid flow, with velocity U that has a sphere of radius R placed in the flow.  There is a obvious axis of symmetry, the ϕ direction, (not to be confused with ϕ[r,θ], the velocity potential) which is oriented parallel to the flow.  The problem will then be solved in the r and θ coordinate directions.  The angle θ = 0, is at the three o'clock position on the sphere and θ=π at nine o'clock (which is the entire range of θ since the angle ϕ would go from 0 to 2 π if necessary -- but not in this problem).

Boundary conditions

Since this is an inviscid flow, the boundary conditions for this flow are quite simple.  As r-> ∞, the velocity must match the free stream flow.  The other boundary condition is that no flow can go through the surface of the sphere.  Our flow is left to right and the angle θ is measured starting at 3 o'clock and increases counter clockwise on the circle.  We can figure out the correct sign for the far away flow by considering that at 9 O'Clock, the flow is purely radial and is coming inward which is negative.  At 3 o'clock the flow is going straight out, which is positive.  Similarly, at 12 o'clock, the flow is purely tangential, in the direction of decreasing θ value.  Thus we use a Cos in front of the radial velocity and a -Sin in front of the tangential velocity.  (I could have drawn a picture, but this is a useful visualization exercise.)    

Thus

For r->∞,
(Subscript[v, r]) = U Cos(θ)
(Subscript[v, θ]) = - U Sin(θ)

For r = R,
(Subscript[v, r]) = 0

Solution using the velocity potential

Governing equation for the flow field in the form of a velocity potential

The governing equation for inviscid flows in three dimensions is Laplace's equation for the potential function, ϕ[r,θ]

Note that this first cell is inactive, but shows the equation, the next one really works

laplaceeq=((((∂((r^2) sin(θ) ((∂ϕ (r, θ))/(∂r))))/(∂r))+((∂(sin(θ) ((∂ϕ (r, θ))/(∂θ))))/(∂θ)))/((r^2) sin (θ)))

laplaceeq = Simplify[1/(r^2*Sin[θ])*
   (D[r^2*Sin[θ]*D[ϕ[r, θ], r], r] + 
     D[Sin[θ]*D[ϕ[r, θ], θ], θ])]

Solution of LaPlace's equation

Linear pde's can often be solved using a "separation of variables" technique. This is done by assuming ϕ[r,θ]  = A[r] G[θ].  As for Re=0 flows, we can often find the form of the angle function from the boundary conditions, in this case, the free stream flow.   

The boundary conditions for the sphere at infinity are (Subscript[v, r]) = U cos[θ], (Subscript[v, θ]) = -U sin[θ].  Recall the definition - ∇ ϕ = v,   thus:  (Subscript[v, r]) =-(( ∂ϕ[r, θ])/∂)r vth = - (1/r) \(\(\(\ \)\(∂ϕ[r, θ]\)\)\/∂θ\) this suggests that the boundary conditions on  ϕ[r,θ] are  
ϕ[r->∞,θ] = -U r Cos[θ],  
at r=R, we have (Subscript[v, r]) = 0, ϕ'[r=R] = 0

This suggests a trial form of phi as

phitrial = A[r] Cos[θ];

We need to start the solution by substituting the trial solution into
the laplace equation

laplacetemp = Simplify[laplaceeq /. 
   {Derivative[a1_, a2_][ϕ][r, θ] :> 
     D[phitrial, {r, a1}, {θ, a2}], ϕ[r, θ] -> phitrial}]

Things are going well, the θ dependence divides out!  Let's get rid of it and rearrange the equation.

laplcesimp = Expand[(laplacetemp/(Cos[θ]/(r^2)))]

This can be solved readily with a solution of the form A[r] = rα because it is an Euler equation.

artemp = 
      A[r] /. DSolve[laplcesimp == 0, A[r], 
            r][[1]]

Here is the far away boundary condition, ϕ[r,θ] --> - r Cos[θ] U (or (Subscript[v, r]) = U Cos[θ])  , @ r-->infinity

bc1 = (artemp /. {C[2] -> 0}) == -r U

The boundary condition on the sphere must be (Subscript[v, r])[r=R] = 0

bc2 = (-(∂_r) artemp /. {r -> R}) == 0

Now solve for the constants

consts = Solve[{bc1, bc2}, {C[1], C[2]}]

phianswer = (artemp /. 
            consts[[1]]) Cos[θ]

We can compute the velocity components  (Subscript[v, r]) =-(( ∂ϕ[r, θ])/∂)r vth = - (1/r) (( ∂ϕ[r, θ])/∂)θ

(Subscript[v, r]) = - D[phianswer,r]

(Subscript[v, θ])  =  - D[phianswer, θ]

Plot of velocity potential for flow past a sphere

Now we can plot the answer in a contour plot.

phifield = 
        ContourPlot[phianswer /. {r -> (Sqrt[(x^2) + (y^2)]), θ -> ArcTan[(y/x)]} /. {R -> 1, U -> 2}, {x, -3, 3}, {y, 0, 
            3}, Contours -> 15, AspectRatio ->  .5];

If you want to see where the sphere is in the picture we can make one.

sphere = 
        ContourPlot[(x^2) + (y^2), {x, -1, 1}, {y, 0, 1}, 
          Contours -> 1, AspectRatio ->  .5];

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.

Velocity potential field for flow past a sphere

Show[phifield, sphere, AspectRatio ->  .5];

We see that potential lines get closer together as the velocity increases.  Further, the potential lines are curved when the velocity gradient is non zero.

Solution using the stream function

Governing equations for the stream function

Now we would like to get the identical answer by using the stream function.  The stream function can always be defined in two-dimensions from the continuity equation and for axisymmetric 3-D flows.  For spherical coordinates, the velocity components are (in terms of ψ)

(Subscript[v, r]) = (1/(((r^2) Sin[θ]))) ∂ψ[r, θ]/∂θ   (Subscript[v, θ]) = -(1/((r Sin[θ]))) ∂ψ[r, θ]/∂r

vrpsi = (∂_θ) ψ[r, θ]/((r^2) Sin[θ]);

vthpsi = -(∂_r) ψ[
                  r, θ]/(r Sin[θ]);

To get  the governing equation for ψ, we use our definition obtained from the continuity equation and substitute into the vorticity components.  Inspection of the components reveals that the (Subscript[Ω, r]) and (Subscript[Ω, θ]) ==0 identically, the only interesting component is (Subscript[Ω, ϕ])  (i. e. coordinate ϕ).  This component is

(Subscript[Ω, ϕ]) = (∂_r)(((r vthpsi))/r) - (∂_θ) (vrpsi/r)

We can simplify this to be

omegafeetemp = 
      Expand[Simplify[-((r (Subscript[Ω, ϕ]))/Csc)[θ]]]

which = 0 if there is no vorticity.  Thus with the restriction of no vorticity, which should not be produced by an inviscid flow, then a simple PDE can be obtained for the stream function.

Note that this is NOT the Laplace equation (It would be for a planar flow!!), even though it looks similar.  Compare to the Laplace equation from above.

Expand[laplaceeq]

Solution of equation by separation of variables

Now we need a functional form for ψ to solve the PDE.  Look at vrpsi above and realize that we need it to have a cos dependence.  Thus C[r] Sin[(θ^2)] looks like a good choice.

psitrial = C[r] Sin[θ]^2;

psieqtemp=omegafeetemp/.{Derivative[a1_, a2_][ψ][r, θ]:>(∂_({r, a1}, {θ, a2}))psitrial,psi[r, θ] -> psitrial}

Things are going well, the theta dependence divides out

\(psieqsimp = Expand[psieqtemp\/\(Sin[θ]\^2\/r\^2\)]\)

This can be solved readily with a solution of the form A[r] = rα because it is an Euler equation.

crtemp = 
      C[r] /. DSolve[psieqsimp == 0, C[r], 
            r][[1]]

Here is the far away boundary condition, we match the r velocity to the free stream , @ r-->∞

bc1psi = ((∂_θ)(((crtemp Sin[θ]^2))/((r^2) Sin[θ])) /. {C[1] -> 0}) == U Cos[θ]

The boundary condition on the sphere must be (Subscript[v, r])[r=R] = 0

bc2psi = ((∂_θ)(((crtemp Sin[θ]^2))/((r^2) Sin[θ])) /. {r -> R}) == 0

Now solve for the constants

consts = Solve[{bc1psi, bc2psi}, {C[1], C[2]}]

Now replace the constants in the answer to get the expression for the streamlines.

psianswer = (crtemp /. 
            consts[[1]]) Sin[θ]^2

Check of potential and stream function answers

Now check the answer from the first method compared to this one.  Here is the radial velocity.

vrpsitest=Simplify[vrpsi/.{Derivative[a1_, a2_][ψ][r, θ]:>(∂_({r, a1}, {θ, a2}))psianswer,psi[r, θ] -> psianswer}]

Compare this to them potential flow answer from above.

Simplify[(Subscript[v, r])]

Now check the tangential velocity

vthpsitest=Simplify[vthpsi/.{Derivative[a1_, a2_][ψ][r, θ]:>(∂_({r, a1}, {θ, a2}))psianswer,psi[r, θ] -> psianswer}]

Compare this to the answer from above.

Simplify[(Subscript[v, θ])]

We are happy to see that they match!!

Plot of streamlines for this flow

psifield = 
        ContourPlot[psianswer /. {r -> (Sqrt[(x^2) + (y^2)]), θ -> ArcTan[(y/x)]} /. {R -> 1, U -> 2}, {x, -3, 3}, {y, 0, 
            3}, Contours -> 15, AspectRatio ->  .5];

Show[psifield, sphere, AspectRatio ->  .5];

Which are the stream lines and they look as we expect.  Again it is important to note that there is no separation on the downstream side of the sphere (which side is it??)

Comparison of streamlines and velocity potential

psifield = 
        ContourPlot[psianswer /. {r -> (Sqrt[(x^2) + (y^2)]), θ -> ArcTan[(y/x)]} /. {R -> 1, U -> 2}, {x, -3, 3}, {y, 0, 
            3}, Contours -> 15, AspectRatio ->  .5, 
          ContourShading -> False];

phifieldx = 
        ContourPlot[phianswer /. {r -> (Sqrt[(x^2) + (y^2)]), θ -> ArcTan[(y/x)]} /. {R -> 1, U -> 2}, {x, -3, 3}, {y, 0, 
            3}, Contours -> 15, AspectRatio ->  .5, 
          ContourShading -> False];

Without shading this looks pretty good.  ψ and ϕ are (trying to be) perpendicular.

Show[sphere, psifield, phifieldx, 
        AspectRatio -> Automatic];

Pressure distribution

Calculation of the pressure field

Now we can calculate the pressure field.  This could be from the Bernoulli equation or from the r component of the N-S equation.  Lets use the latter.

dpdr = -(vrpsi*D[vrpsi, r] + (vthpsi*D[vrpsi, θ])/r - 
    	vthpsi^2/r)

dpdrtemp = Simplify[dpdr /. 
    {Derivative[a1_, a2_][ψ][r, θ] :> 
      D[psianswer, {r, a1}, {θ, a2}], 
     ψ[r, th] -> psianswer}]

If we form the definite integral for pressure as a function of r, we would have p(r->∞) - p(r) = (∫_r^∞)∂(p/∂)r

pofrtemp = -(∫_r^rwayway) 
              dpdrtemp \[DifferentialD]r +  pfaraway

It is clear that these die off in distance.  Let's make the far away pressure = 0, i.e., rwayway ==> really far away.  If this is the case, then the second and 4th terms are 0.   The pressure as a function of r and theta is:

pofr = 
      Simplify[Limit[(pofrtemp /. pfaraway -> 0), 
          rwayway -> ∞]]

It is interesting to see what the pressure on the surface of the sphere is.

psurface=Simplify[pofr/.r->R]

Now what is the maximum pressure?

pstag=psurface/.θ->0

Plot of the pressure field

pfield = 
        ContourPlot[pofr /. {r -> (Sqrt[(x^2) + (y^2)]), θ -> 
                  ArcTan[(y/x)]} /. {R -> 1, U -> 2}, {x, -5, 
            5}, {y, 0, 5}, Contours -> 10, 
          AspectRatio ->  .5];

Show[pfield, sphere, AspectRatio ->  .5];

Other issues

Flow separation

The pressure field is symmetric.  The highest pressure is in front of and behind the flow in the stagnation point region.  The lowest pressure is at the top where the velocity is fastest.  From the pressure field, we can see that separation is a likely possibility for a real fluid.  Consider a left to right flow.  For x<0, the value of the pressure is always decreasing, this is a "favorable" pressure gradient.  If you were a fluid particle, you would feel more push behind than in front.  For x>0, the pressure is increasing.  The fluid encounters an "adverse" pressure gradient and the fluid is feeling more push in front than behind and this slowing down.  If we remember that there is really drag caused by the solid surface, a fluid particle will not have enough momentum to get from the π/2 point back down to θ=0 (as it is expected to from the no drag case).  Thus the fluid has no choice but to go off at some angle up to the right.

Drag on the sphere

There are two kinds of drag that are defined.  One is "skin drag" which is the integral of the tangential force on the surface of the sphere. This is defined to be 0 for inviscid fluids.  The "form drag" is the integral of the total normal stress over the entire sphere.  Because the pressure is the only normal stress for an inviscid fluid, we can guess that the form drag is likewise 0 because the picture is clearly symmetric.

Let's calculate the form drag anyway.   See Denn page 253 for the low Re case.

drag = (∫_0^(2 π))((∫_0^π) (R^2) Sin[θ] (pofr /. 
                    r -> 
                      R) Cos[θ] \[DifferentialD]θ) \[DifferentialD]ϕ