Boundary-layer flow over a flat plate and wedge

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:  11/30/17

Summary

This notebook examines the boundary-layer flow over a flat plate and a wedge.  For the flat plate, a similarity variable is used to reduce the continuity and momentum equations to a single nonlinear ODE. (A slightly more complex transformation and corresponding ODE exist for the flow over a wedge.)  For both cases, the solution is given using a numerical scheme and several important physical aspects of boundary layers are elucidated.

References

The reference for this NoteBook is:
M. M. Denn
Process Fluid Mechanics
Prentice Hall
1978

R. L. Panton
Incompressible Flow
Wiley
1986

Homework problem for this note book

1.  How does the boundary layer thickness change with Reynolds number ??
2.  How does the wall shear stress change with Reynolds number ??
3.  What is the primary physical characteristic of a boundary layers ??
4.  Run the code to find is the value of η for which the velocity is 95% of the free stream. 
5.  Run the code to find the value of the stress at the wall for a flat plate.  Can you get more significant figures?  Check elsewhere to see if my value of f’’’(0) == .332058 matches other calculations.
6.  Run the code to show how the boundary-layer thickness and shear stress change when the  imposed pressure gradient changes (i.e. the wedge problem)  Find the angle at which separation occurs.  Look up the solutions for flow around a cylinder.  How does the angle for separation for the wedge compare to the angle of separation for flow around a cylinder (from “Boundary-Layer Theory, by H. Schlichting, p161, this occurs at θ = 108.8 as measured from the side that the flow impinges on the cylinder)  Note that these angles are measured differently so make a sensible comparison.
7.  The Reynolds numbers for cars, planes and boats is very large.   What is the boundary layer thickness in real variables  (say cm) 25 cm from the leading edge of
a.  A WWII Spitfire airplane
b.  A 787 landing
c.  An “8’s” racing shell
d.  Your first new car cruising to your new job at 55 mph!
8.  For different values of the wedge angle, find wall shear stress and compare the exact result to just using τ=(U/δ)
9.  Explain the relationship between the variable for the suddenly started plate: y/(Sqrt[ν t]),  and η =(y/(Sqrt[((x ν)/U)]))for steady flow past a flat plate.  Note that ν is “nu” even if the font is not clear.

Flow past a flat plate

Do background work to set up the problem

Existence of a similarity variable.

A good reference for this section is L. G. Leal  (1992) Laminar Flow and Convective Transport Processes,  Butterworth.)  

If you encounter a fluid flow problem in a rectangular channel, you expect the characteristic length to be the channel height.  For a circular pipe, it will be the diameter or radius.  For flow past a sphere, the radius is a reasonable length.  Thus we always expect to find a characteristic length, that is contained in the physics of the flow.  However, for a boundary-layer flow, where the boundary-layer is much smaller than the scale of the container, or a boundary-layer on the outside of an object where again the layer is much (say ~1/100) smaller than the device, what now is the relevant length scale?  Would nature leave us with no characteristic length?  

The answer is that there is a natural length scale.  It is just not a geometric length scale.  In the absence of a reasonable geometric length scale we look to define a "similarity" variable.  The implication is that the velocity or temperature profiles are "self-similar", and with this proper scaling can be collapsed onto the same curve.  The length scale for the normal direction within a boundary - layer is constructed from the physical properties of the flow and any other length or time variables that come to play in the problem.

For example, for the "suddenly started" plate problem, where a long and wide flat plate is started in motion, the boundary-layer thickness scales as y/(Sqrt[ν t]), where t is time and ν is the kinematic viscosity. (The boundary layer gets thicker in time as (Sqrt[ν t]).) 

To wrap this up, you always want to investigate the scaling of your problems.  If there is no geometric length scale, (i.e., problem has a infinite spatial domain), you expect the existence of a similarity variable which is a combination of the other variables and physical properties.  

Use of a similarity variable can lead to transformations that turn the PDE's into ODE's.  This can work for both linear and nonlinear equations.

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

Boundary layers are regions where two competing transport effects are about the same order of magnitude.It is expected that the gradients are also "large" in this region.  The scaling arguments on pp 189-191 in Truskey et al.  may also help us understand the simplifications that reduced the momentum equation to:

(vx ((∂(Subscript[v, x]))/(∂x))+vy ((∂(Subscript[v, x]))/(∂y)))=-(1/ρ)((∂P)/(∂x))+ν(((∂^2)(Subscript[v, x]))/(∂y2))

If the plate is flat there is no change in pressure because of a change in the fluid local velocity,say for example if the flow were converging or diverging.  You might have a slight pressuredropin a closed system,but this can usually be neglected.Therefore,((∂P)/(∂x)) =0.

The continuity equation is just what it always is.

conteq = D[u[x,y],x] +D[v[x,y],y]

The momentum equation is

umomeq = u[x,y] D[u[x,y],x] + v[x,y] D[u[x,y],y] - ν   D[u[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.

Note that continuity equation is satisfied automatically by a stream function formulation,
 u  = ((∂ψ[x,y])/(∂y)) ,  v  =- ((∂ψ[x,y])/(∂x)).

conteq/.{u[x,y]->D[ψ[x,y],y],v[x,y]->-D[ψ[x,y],x],Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))D[ψ[x,y],y],
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))D[-ψ[x,y],x]}

Transformation of the momentum equation in terms of the stream function

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

psibleq=umomeq/.{u[x,y]->D[ψ[x,y],y],v[x,y]->-D[ψ[x,y],x],Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))D[ψ[x,y],y],
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))D[-ψ[x,y],x]}

We won't use this equation because the substitutions are harder.

Similarity variable formulation

Because the only length scale is the one measured along the plate and because we expect that the velocity profiles will be self similar (the same shape scaled by just the local boundary-layer thickness), we look for similarity variables.

Start by defining a similarity variable that will be dimensionless.  We will just choose the correct one (see Panton or Denn), there are formal procedures for finding these.  (See Leal's book for how to do this.)  We will also be changing the dimensional equation into a dimensionless ODE.  When we do this, there will not be a parameter in this resulting equation because the length scale for y is the boundary-layer thickness (which scales as (Re^(-1/2))).  Thus the introduction of a second length scale which contains the Reynolds number removes the Reynolds number in the equation.

eta=y /(Sqrt[ ν x/U])

η ≡(y/(Sqrt[((x ν)/U)]))

Likewise for the stream function,

psi= f[η[x,y]] (Sqrt[ ν x U])

We use the Greek, in the equations, but for substitution purposes, it is convenient to have a different name, eta and psi.

Define the velocities in terms of the stream function and similarity variable

We need to use the chain rule to find all of the derivative substitutions and to find (Subscript[v, x]) and (Subscript[v, y]). 

First find (Subscript[v, x])(x,y),  (Subscript[v, x])  = ((∂ψ[x,y])/(∂y)) =   ((∂ψ)/(∂η))((∂η)/(∂y))

We want the velocity in terms of our definition of the stream function.  Mathematica knows the chain rule.

eu1=  D[psi,y]

We can also get (Subscript[v, x]) by using the chain rule explicitly and fortunately they are the same...

eu1= D[psi,η[x,y]] D[η[x,y],y]

We can simplify this by using our similarity variable to give (Subscript[v, x])/U = f'[η], where U is the free stream velocity.

eu2 =PowerExpand[ eu1/.D[η[x,y],y]->D[eta,y]]

Note again that U shows up because the f(η) is dimensionless and the flow direction velocity is scaled with U.

Now get the y derivative of η and substitute.

Now we need v.
 (Subscript[v, y])  =- ((∂ψ[x,y])/(∂x))  =  - ((∂ψ)/(∂η))((∂η)/(∂x))

vee1=-D[psi,x]

vee2=vee1/.D[η[x,y],x]->D[eta,x]

For further substitutions below we it is easiest to use vee2, however we can simplify this some more to make it look nice,

vee3=PowerExpand[vee2/.η[x,y] ->η]

We now have (Subscript[v, y])/(Sqrt[x/Uν])=1/2 (η f'[η]- f[η])

vee4=PowerExpand[vee3/.y->η (Sqrt[ ν x/U])]

Now derive the nonlinear ODE that describes the flow field.

Here is the boundary-layer momentum equation

umomeq

Now we make the substitutions for the velocities,

ueq1=umomeq/.{u[x,y]->eu2,v[x,y]->vee2,Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))eu2,
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))vee2}

This step replaces derivatives of η.

ueq2= ueq1/.{
D[η[x,y],x]->D[eta,x],D[η[x,y],y]->D[eta,y],
D[η[x,y],{x,2}]->D[eta,{x,2}],D[η[x,y],{y,2}]->D[eta,{y,2}]}

Now that we are done taking derivatives we can make the η(x,y) into just η.

ueq3=ueq2/.η[x,y]->η

We now clean it all up and remove common factors,

ufinal=ExpandAll[-PowerExpand[ueq3 /U^2x]]

Return to conclusions

The final boundary-Layer ODE is

ufinal

Boundary conditions

The boundary conditions are that there is no slip on the surface, no flow through the surface and that the tangential velocity matches the free stream far away. 

The are in terms of the newest notation:

u(x,y=0)=0, η = 0, f'(η= 0)=0,
v(x,y=0)=0, η = 0, f(η= 0=0,
y->∞, u(x,y->∞)=U, f'(η=∞) = 1,  
x=0, u(x,y->∞)=U, f'(η=∞) = 1,

The primes (') denote derivatives with respect to the similarity variable.  This nonlinear ode does not have a known analytical solution.  However, it can be easily solved numerically.  
		
Numerical solutions for ode's are often done by creating a system of first order odes.  This is done by defining
		
		y1 = f,
		y2 = f'
		y3 = f'',
		
These have to be related. 
		
		Thus we have
		y1' = y2 (by definition)
		y2' = y3 (by definition)
		y3' = -1/2*y1*y3   (which is from the original ode.)
		
In the first part I solve the flat plate problem separately using a crude numerical scheme that usually converges, although not real fast.  It is a shooting method.  I have picked y =15 as the end of the integration.  The NDSolve routine uses a Runge-Kutta method with built in step size adjustment.  

Note that “x” in this calculation is the variable η in our analysis.

fppinit=1;
eps=1;

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

Here is the tangential velocity

Plot[y2[x]/. zz[[1]],{x,0,10},PlotRange->All,AxesLabel->{"η","\!\(\*SubscriptBox[\(v\), \(x\)]\)/Uo"}]

Here is the dimensionless normal velocity.   Note also that it does not go to 0 at the top of the layer.  Fluid is being expelled as the layer thickens!

Plot[(1/2) (x y2[x]-y1[x])/. zz[[1]],{x,0,10},PlotRange->All,AxesLabel->{"η","\!\(\*SubscriptBox[\(v\), \(y\)]\)/\!\(\*SqrtBox[\(νUo/x\)]\)"}]

The second derivative, f'' is related to the stress on the plate.

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

Three important results are obtained.  First, the graphs show that the effective boundary layer thickness is about 6 in η.  Second, the vertical velocity is not 0 at the top of the boundary layer -- fluid is being expelled.  It has to be, the fluid inside the boundary layer is being slowed down.
Third, the wall shear stress and other transport properties are related to the value of f''[0] = .332058.

Here are tables of values

Table[y1[x]/.zz[[1]],{x,0,8.8,.2}]

U* y2[η] is the flow direction velocity

Table[y2[x]/.zz[[1]],{x,0,8.8,.2}]

y3 is related to the shear stress.

Table[y3[x]/.zz[[1]],{x,0,8.8,.2}]

Flow past a wedge

Governing equations and problem set up

From Denn the flow past a wedge is given by
U(x) = A xm  , m=beta/(2 pi-beta)

Recall that the pressure gradient, ((∂p)/(∂x)), will change as -U(x) U'(x). This gives

 (vx ((∂(Subscript[v, x]))/(∂x)) + vy ((∂(Subscript[v, x]))/(∂y))) = U(x) U'(x) + ν (((∂^2)(Subscript[v, x]))/(∂y2))
 
 
******  Thus, the flow past a wedge is a model problem for telling how the boundary layer will change as the pressure gradient changes. *****

Denn tells us that the equation can be reduced to a nonlinear ode again taking advantage of similarity of solution profiles.  This gives

		f''' + (m+1)/2 f f'' + m (1-f'2)  = 0.
		
		If m=0, the equation is the same as for the flat plate.  


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

y1' = y2
y2' = y3
y3' = - (m+1)/2 y1 y3 + m*(1-y22)   

This time I use a real shooting method, with a pseudo Newton-Raphson iteration.  It works for most cases.  If the angle is too negative, the stress will go through 0, this means that there is no layer and thus it is not surprising that the answer blows up.  Of course if we look at the physical situation, we might be surprised that we can get any solution for negative β.  The trend is correct and it is interesting to examine the case of a diverging flow.

Numerical solution

Note that you will want to run this several times changing the value of β over the greatest possible range.  If β=0, the result is the same as a flat plate.

fppinit=1;
eps=1;
beta=-.621;
m=(beta/(2 π-beta));
youter=10;

angleindegrees=N[((beta 180)/π)]

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]==-(1/2) (m+1) y1[x] y3[x]-m (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]>.000001,zz=NDSolve[{y1'[x]==y2[x],y2'[x]==y3[x],y3'[x]==-(1/2) (m+1) y1[x] y3[x]-m (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->{"zeta","\!\(\*SubscriptBox[\(v\), \(x\)]\)/Uo"}]

Again note that this is magnitude of the normal velocity after it has been rescaled to be of order 1.  It is really smaller by the relative thickness of the boundary layer to the distance variable -- as required by the continuity equation.

Plot[(1/2) (x y2[x]-y1[x])/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"zeta","\!\(\*SubscriptBox[\(v\), \(y\)]\)/\!\(\*SqrtBox[\(νUo/x\)]\)"}]

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

Plots of the results at different angles

Here I save some old results, which I must do after a run. (You need to do these yourself, they will not just run !!).  Just make up a name for the result that you wish to save and then equate this to "zz" as I show below.  Your newly chosen name will store all of the useful results of the calculation.

For -.622I did

zzmp622=zz;

For -.6 I did

zzmp6=zz;

For -.22 I did

zzmp22=zz;

For 0

zz0=zz;

-.3

zzmp3=zz;

.4

zzp4=zz;

1

zz1=zz;

Here I plot some of them together.  The middle one is the flow plate.  The two to the left of it are for a converging flow (more stress, thinner layer), the two to the right are for diverging flow (thicker layer, less stress).

Plot[{y2[x]/. zzmp622[[1]],y2[x]/. zzmp6[[1]],y2[x]/. zzmp3[[1]],y2[x]/. zz0[[1]],y2[x]/. zzp4[[1]],y2[x]/. zz1[[1]]},{x,0,youter},PlotRange->All,AxesLabel->{"η","f'"},PlotStyle->{Dashing[{0.02`,0.02`}],Dashing[{0.05`,0.03`}],Dashing[{1,0}],Dashing[{0.03`,0.01`}],Dashing[{0.02`,0.04`,0.06`,0.04`}]},PlotLegends->{"β=-.622","β=-.6","β=-.3","β=0","β=.4","β=1"}]

Plot[{y3[x]/. zzmp622[[1]],y3[x]/. zzmp6[[1]],y3[x]/. zzmp3[[1]],y3[x]/. zz0[[1]],y3[x]/. zzp4[[1]],y3[x]/. zz1[[1]]},{x,0,youter},PlotRange->All,AxesLabel->{"η","f''"},PlotStyle->{Dashing[{0.02`,0.02`}],Dashing[{0.05`,0.03`}],Dashing[{1,0}],Dashing[{0  .03`,0.01`}],Dashing[{0.02`,0.04`,0.06`,0.04`}]},PlotLegends->{"β=-.622","β=-.6","β=-.3","β=0","β=.4","β=1"}]

The important result is that in a converging region, β>0, the stress increases and the layer thins.  In a diverging region, β<0, the stress decreases, to 0, and the layer thickens.   This tells us the qualitative behavior for flow past any shape body.  Further it provides insight into boundary layer separation which occurs when the tangential velocity gradient at the wall equals zero indicating that a backflow will start to occur.