Laminar boundary-layer flow past a flat plate

>

Consider a fluid stream with velocity and kinematic viscosity in which a thin plate is inserted parallel with the fluid flow. Determining the velocity of the fluid in the region close to the plate is the Blasius problem [6, p. 233]. Assuming the flow is steady, incompressible, and Newtonian, the plate is infinitely wide, and neglecting buoyancy, the equations of motion and continuity are:

> restart;

> alias( U=u(x,y), V=v(x,y) ):

> PDE:={ U*diff(U,x)+V*diff(U,y)-nu*diff(U,y$2)=0,

> diff(U,x)+diff(V,y)=0 };

where and are the - and -components of the fluid velocity. The boundary conditions consist of the ``no-slip'' conditions on the boundary of the plate: = =0, and the free stream-merge condition , for all >0.

>

A similarity transformation can be used to reduce this parabolic system of PDEs to a single ODE. This can be done by choosing the dimensionless similarity variable to be:

> simsubs:=eta(x,y)=y*sqrt(u[0]/nu/x/2);

The corresponding nondimensional stream function for the flow is:

> stream:=psi(x,y)=sqrt(2*nu*x*u[0])*f(eta(x,y));

Thus, the velocities can be expressed as:

> Usubs:=U=diff(rhs(stream),y);

> Vsubs:=V=-diff(rhs(stream),x);

Substituting the stream function representations of the velocities into the PDEs is tedious to complete by hand. Fortunately, this is exactly one of Maple's strengths:

> ODE:=simplify(subs(Usubs,Vsubs,simsubs,PDE));

The trivial fulfillment of the continuity equation is evident in this result. However, the first equation is not so readily identified. To simplify this further, note that each argument to , and its derivatives, is simply . To force this simplification,

> simsubs2:=solve(subs(eta(x,y)=eta,simsubs),{y});

> ODE:=simplify(subs(simsubs2,ODE),symbolic);

It is now easy to identify the Blasius equation

, >0.

The conversion of the boundary conditions can be done by inspection. The resulting conditions are = = and .

>

The ``boundary condition'' at presents a problem. However, a simple asymptotic analysis shows that solutions at =10 are safely in the far-field. The problem now takes the form of a third-order two-point boundary value problem for on [0,10]. The reformulation as a first-order system is straightforward; let and , then

,

, (4)

,

where is the control parameter. The objective is to find so that the solution to (4) satisfies the boundary condition . Since there is only one boundary condition at we have and . The sensitivity equations that are obtained from (4), together with their accompanying initial values, are:

,

, (5)

,

where, for example, . Now, assume that an approximate solution to (4) and (5) has been computed with , i.e. , . Then, assuming is not sufficiently small, the next iteration will be made with

.

>

While it is important to understand how shoot obtains its approximations, it is also a tremendous advantage to use Maple to both compute automatically, and symbolically, the sensitivity equations and iteratively solve the combined system of IVPs and take one step of the Newton-Raphson method until the boundary conditions are satisfied within a prescribed tolerance.

>

Prior to preparing this problem for solution by shoot , it is necessary to load the definition of this procedure:

> with(Shoot):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the assigned name shoot now has a global binding

Now, define the dependent variables for the first-order system:

> FNS:={ f(eta), g(eta), h(eta) }:

and the differential equations which they satisfy:

> ODE:={ diff(f(eta),eta)=g(eta),

> diff(g(eta),eta)=h(eta),

> diff(h(eta),eta)=-f(eta)*h(eta) }:

The initial condition for is unknown, and the boundary condition at is the control for our iterations:

> IC:={ f(0)=0, g(0)=0, h(0)=beta }:

> BC:=g(10)=1:

We take the first shot with , and iterate until the boundary condition is satisfied to six decimal places:

> S:=shoot( ODE, IC, BC, FNS, beta=0,

> abserr=Float(5,-7), output=listprocedure, method=taylorseries ):

This solution can be analyzed using any of the standard Maple tools for manipulating the numerical solution of a differential equation. A common place to begin is a plot of the solution, including labels to distinguish the different curves:

> with(plots):

> C:=odeplot( S, [ [eta,f(eta)],

> [eta,g(eta)], [eta,h(eta)] ],0..4):

> subs( [g=`f '`, h=`f ''`],

> subs( [seq( op(0,lhs(F))=rhs(F(3.5)),

> F=subsop(1=NULL,S) )] ) ):

> L:=textplot(

> subs( eta=3.5,

> map(F->[eta,rhs(F)+1/10,lhs(F)],%)),

> font=[TIMES,ITALIC,18], align=ABOVE ):

> display( {C,L}, labels=[`eta`,``],

> title=`Blasius solution for flat-plate `||

> `boundary layer` );

Warning, the name changecoords has been redefined

With increasing distance from the leading edge of the plate in the downstream direction the thickness of the retarded boundary layer increases continuously as increasing quantities of fluid become affected. In the boundary layer the velocity of the fluid increases from zero at the wall (no slip) to its full value which corresponds to external frictionless flow. The velocity reaches 99% of its bulk value when :

> fp0:=subs( S, g(eta) ):

> fp := proc(x)

> if not type(evalf(x),`numeric`) then

> 'procname'(x);

> else

> fp0(x);

> end if;

> end proc:

> eta[`99%`]:=fsolve( fp(eta)=0.99, eta=3..4 );

Note that the computation of has changed significantly for Maple 6. The solution returned by shoot is obtained with method=taylorseries and the function used in fsolve must be written to avoid evaluation in cases when the argument is not numeric. But, in the end, the result is equivalent to the one reported in the original paper.

>

The corresponding boundary layer thickness is:

> delta[`99%`]:=subs(eta=eta[`99%`],

> combine(rhs(op(simsubs2))) );

As is now evident, the thickness of the boundary layer decreases with decreasing viscosity.

>

The shear stress, , is:

> tau:=mu*simplify(diff(subs(simsubs,rhs(Usubs)),y));

Note that a large velocity gradient across the flow creates considerable shear stress in the boundary layer. The wall shear stress is:

> tau[w]:=subs( (D@@2)(f)(0)=subs(S,h(eta))(0), subs( y=0, tau ) );

The analysis can be continued to obtain an understanding of parameters such as the displacement thickness and drag coefficient.

>