Inviscid flow past a stationary cylinder

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

Version:  4/11/16

Summary

Here we consider the problem of a uniform planar flow encountering a stationary circle.  That is, a constant velocity flow past an infinitely long circular cylinder.

Governing equation for the flow field

We discussed in class that the governing equation for inviscid flows in a plane dimensions is Laplace's equation for the stream function, ψ[r,θ]

laplaceeq = 1/r*D[r  D[ψ[r,θ],r],r]+1/r^2 D[ψ[r,θ],{θ,2}]

Linear pde's can often be solved using a "separation of variables" technique. This is done by assuming ψ[r,θ] = f[r] g[θ].  As for the Re=0 flows, we can find the form of the angle function from the boundary conditions.   

The boundary conditions for the sphere at infinity are vr = U cos[th], vθ = -U sin[th].  Recall the definition - grad phi = v,   thus:  vr =- D[phi[r,th],r], vth = - 1/r D[phi[r,th],th] this suggests that the boundary conditions on phi[r,th] are phi[r-->inf] = -U r Cos[th],  
at r=R, we have vr = 0, phi'[r=R] = 0

This suggests a trial form of phi as

ψtrial = f[r] Sin[th] ;

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

laplacetemp= laplaceeq/.{Derivative[a1_,a2_][ψ][r,θ]:>
		 D[ψtrial,{r,a1},{th,a2}],
		 ψ[r,th]->ψtrial}

Things are going well, the theta dependence divides out.

laplcesimp = Expand[laplacetemp/(  Sin[th]/r^2 )]

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

frtemp = f[r]/.(DSolve[{laplcesimp==0,f[R]==0},f[r],r] )[[1]]

frtemp1 =Expand[frtemp]

We see that we want the first term to match U r

frtemp2 = Solve[frtemp1 [[1]] == U r, C[1]]

So the answer for f[r] is

frfinal = frtemp/.frtemp2[[1]]

ψanswer = frfinal Sin[θ]

Now we can plot the answer in a contour plot.

streamlines1=ContourPlot[(ψanswer/.{r->Sqrt[x^2+y^2],
θ->ArcTan[y/x]}/.{R->1,U->2}),{x,0,5},{y,0,5},RegionFunction->Function[{x,y},x^2+y^2>1],
Contours->20,BoundaryStyle->Red, AspectRatio->.5]

Note that the contours smoothly change as we approach the sphere.  The total flow between each streamline is constant.

streamlines2=ContourPlot[(ψanswer/.{r->Sqrt[x^2+y^2],
θ->ArcTan[y/x]+π}/.{R->1,U->2}),{x,-5,0},{y,0,5},RegionFunction->Function[{x,y},x^2+y^2>1],
Contours->20,BoundaryStyle->Red, AspectRatio->.5]

Show[streamlines1,streamlines2,PlotRange->All]

It is easy to see that at θ=π/2, the velocity has increased compared to its value at θ=0 or π/2.

The velocity components are:

vrr = 1/r D[ψanswer,θ]

vθθ =Simplify[- D[ψanswer,r]]

Find the pressure field at r =R, P(r->infinity) ==>0

The Bernoulli equation gives

pθR =Simplify[ U^2/2- vθθ^2/2]

pressure1=ContourPlot[(pθR/.{r->Sqrt[x^2+y^2],
θ->ArcTan[y/x]}/.{R->1,U->2}),{x,0,5},{y,0,5},RegionFunction->Function[{x,y},x^2+y^2>1],BoundaryStyle->Red, AspectRatio->.5,Contours->10]

pressure2=ContourPlot[(pθR/.{r->Sqrt[x^2+y^2],
θ->ArcTan[y/x]+π}/.{R->1,U->2}),{x,-5,0},{y,0,5},RegionFunction->Function[{x,y},x^2+y^2>1],BoundaryStyle->Red, AspectRatio->.5,Contours->10]

Show[pressure1,pressure2,PlotRange->All]

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 Pi/2 point back down to th=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.

drag=2 R Integrate[ prθ   Cos[θ]  ,{θ,0, π}]

As expected, the drag is 0.