Stagnation flow

This notebook has been written in Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular 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/7/14

Summary

This notebook examines a stagation flow at high Reynolds number.  The other region is inviscid and in the region is a boundary-layer flow.  Because the outer velocity in the x direction is accelerating, or because the pressure in the outer region is decelerating, the boundary layer remains of constant thickness.

References

The references for this NoteBook are:
F. M. White
Viscous Fluid Flow
McGraw-Hill, 1974

G. K. Batchelor
An Introduction to Fluid Dynamics
Cambridge, 1967

Flow past a flat plate

Start with the vocticity transport equation

This “special” case of the general flow past a solid body, for which a similarity variable can be found, is that y and x are related because y = const for all x.

The  mass and momentum  equations for a forced boundary - layer flow past a flat plate

The vorticity transport equation is

(vx ((∂ω)/(∂x))+vy ((∂ω)/(∂y)))= ν((((∂^2)ω)/(∂x2))+(((∂^2)ω)/(∂y2)))

We could use the x and y Navier-Stokes equations, but this is more fun.

The vorticity equation is

vorticityeq = u[x,y] D[ω[x,y],x] + v[x,y] D[ω[x,y],y] - ν   (  D[ω[x,y],{x,2}]+D[ω[x,y],{y,2}])

Stream function formulation to satisfy continuity

We want to turn these two coupled PDE's into a single ODE.  It will still be nonlinear, but it is much more easy to solve that the original 2 equations.

We will define a stream function from the continuity equation, define a similarity variable involving x and y.  Then transform the equations.

How does the vorticity come from the stream function

-D[u[x,y],y] + D[v[x,y],x]

%/.{u[x,y]->B x f'[y],v[x,y]->-B f[y],Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))(B x f'[y]),
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))(-B f[y]),Derivative[a1_,a2_][ω][x,y]:>  (∂_({x,a1},{y,a2}))(-B x f''[y])}

As mentioned before, the stream function ψ, can represent both the x and y directions velocities, this gives for the momentum equations,

vortbleq=vorticityeq /.{u[x,y]->B x f'[y],v[x,y]->-B f[y],Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))(B x f'[y]),
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))(-B f[y]),Derivative[a1_,a2_][ω][x,y]:>  (∂_({x,a1},{y,a2}))(-B x f''[y])}

Expand[%/B/x/ν]

This equation in terms of exact differentials is

B/ν D[f[y] f''[y],y]-B/ν D[ f'[y]^2,y] + f''''[y]

Which matches, as a check.

Integrate [ B/ν D[f[y] f''[y],y]-B/ν D[ f'[y]^2,y] + f''''[y],y  ]

Which we know must equal a constant.

For f(infinity) = 1 means the constant C = - B/ν, which gives...

f'''[y] + B/ν f[y] f''[y] + B/ν (1- f'[y]^2)

Now nondimensionalize this to get rid of B and ν

%/. {Derivative[a1_][f][y]:>  (∂_({η,a1}))(F[η]/Sqrt[B/ν])Sqrt[B/ν]^a1,f[y]->F[η]/Sqrt[B/ν]}

ExpandAll[% ν/B]

Stagnation flow calculations

Governing equations and problem set up

The equation is

		f''' +  f f'' +  (1-f'2)  = 0.
		

We again need to make a system of first order ODE's .  These are

y1' = y2
y2' = y3
y3' = - y1 y3 + β*(1-y22)   

This time I use a real shooting method, with a pseudo Newton-Raphson iteration.  It works for most cases.

Numerical solution

fppinit=1;
eps=1;
youter=4;

need to run one time to get a value for the outer f'[infinity]

zz=NDSolve[{y1'[x]==y2[x],y2'[x]==y3[x],y3'[x]==- y1[x] y3[x]-  (1-((y2[x])^2)),y1[0]==0,y2[0]==0,y3[0]==fppinit},{y1[x],y2[x],y3[x]},{x,0,youter}];
xz=(zz/.x->youter)[[1,2,2]];
eps=1-xz;
Print["U = ",xz,"    error=  ",eps,"   f''[0]=  ",fppinit];

fppinitold=fppinit;
fppinit=fppinit+.01;
xzold=xz;

While[Abs[eps]>.00001,zz=NDSolve[{y1'[x]==y2[x],y2'[x]==y3[x],y3'[x]==- y1[x] y3[x]-β (1-((y2[x])^2)),y1[0]==0,y2[0]==0,y3[0]==fppinit},{y1[x],y2[x],y3[x]},{x,0,youter}];xz=(zz/.x->youter)[[1,2,2]];eps=1-xz;xcorrect=(((1-xz) (fppinit-fppinitold))/(xz-xzold));fppinitold=fppinit;xzold=xz;fppinit=fppinitold+xcorrect;Print["U = ",xz," xzold= ",xzold,"    error=  ",eps,"   f''[0]=  ",fppinit," xc= ",xcorrect];xzold=xz;]

Plot[y2[x]/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{η},PlotLegends->{"f'"},PlotStyle->Red]

So the thickness is about 2.4.  
Here is the normal velocity which is lot different than than a flat plate

Plot[ (y1[x])/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"f"},PlotStyle->Green]

Plot[y3[x]/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"f''"},PlotStyle->Blue]

Show[%,%%,%%%]

Plot[(-1/2 y1[x]^2- y2[x])/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"pressure"},PlotStyle->Orange,PlotLabel->"Pressure drop away from stagnation point"]

Plot[(-1/2 y1[x]^2- y2[x])/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"inner pressure"},PlotStyle->Orange,PlotLabel->"Pressure drop away from stagnation point"]

Plot[-1/2 y^2,{y,7,0},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"outer pressure"},PlotStyle->Green,PlotLabel->"Pressure drop away from stagnation point"]

Show[%,%%]

Plot[ (y1[x])/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"f"},PlotStyle->Green]

Plot[ y/. zz[[1]],{y,7,0},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"outer v"},PlotStyle->Purple]

Show[%,%%]

Plot[ y-.65/. zz[[1]],{y,7,0},PlotRange->All,AxesLabel->{"η"},PlotLegends->{"outer v"},PlotStyle->Purple]

Show[%,%145]

So the effect of the wall is to displace the outer region by 0.65 in y.