BOUNDARY LAYER CODES FOR STUDENTS
William J. Devenport
Joseph A. Schetz
Aerospace and Ocean Engineering Department
Here we present simple JAVA computer codes that are intended for student use solving the homework problems in a text such as "Boundary Layer Analysis" by J. A. Schetz (Prentice Hall, 1993) and other similar problems on a PC or Work Station. These codes are specifically not intended as general purpose codes for use by working professionals in the field. The goal has been to keep the formulation, logic and programming as simple as possible so that the student can easily grasp the flow of the calculations. Thus, primitive variables (u,v); (x,y) are employed with no transformations. Last, only limited physical property information is included. The codes have all been compiled and tested thoroughly.
These simple codes are meant only to relieve the student of the burden of writing and debugging numerous codes during a semester or two course. The hope is to thereby leave sufficient time and energy for the working of actual boundary layer problems of reasonable complexity. The student is not relieved of the burden to think. One should always estimate the answer the computer solution is expected to yield. For example, one could use a simple integral solution to estimate the level of the skin friction to be expected. The student must not fall into the trap of assuming that because the computer produces nice, neat output it must be correct. The student is the analyst, not the computer, and he or she is responsible for producing the correct answer.
There are eight separate codes for the various methods covered: 1) Thwaites-Walz incompressible, laminar integral method, 2) incompressible, laminar boundary layers by an explicit numerical method, 3) incompressible, laminar boundary layers by an implicit numerical method, 4) compressible, laminar boundary layers by an implicit numerical method, 5) Moses, incompressible turbulent integral method, 6) incompressible turbulent boundary layers by an implicit numerical method, 7) incompressible turbulent boundary layers by an implicit numerical method with a stretched grid, and 8) turbulent jets and wakes by an implicit numerical method. A brief description of each code precedes the listings. Also, a default input for a typical example and results for that case are presented and discussed. Detailed development of the methods can be found in the referenced text. Earlier Fortran versions of these codes are also in that reference.
For any boundary layer problem, one must first specify the fluid through density,r and viscosity, m, or perhaps just kinematic viscosity, n. In variable-density, variable-property cases, an equation of state and property variation information must also be given. The next information needed would be the freestream velocity, U¥ , and the inviscid edge velocity distribution, Ue(x), (and sometimes the edge velocity gradient dUe/dx). For high-speed flows, Mach numbers, M¥ and Me(x), are needed. With heat transfer or compressible flow, the wall temperature, TW(x), or the wall heat transfer distribution and T¥ and Te(x) are required. Mass transfer cases are not included in the codes given here. Next is the matter of initial conditions. The integral methods only require initial values of quantities like q (xi) and Cf (xi). The differential methods need complete, initial profiles u(xi,y), v(xi,y). For turbulent cases, one must select the turbulence model from those available in the code, e.g. mixing length model, eddy viscosity model or TKE model with the Prandtl Energy Method. The TKE method requires an initial profile and boundary conditions for K. There are also numerous so-called modeling constants such as k whose values the user may wish to change from the values that come in the code. Also, a transition location must be specified.
Finally, the computational region must be selected. With an integral method, this is simply the length of interest. The differential methods require the height and the length of the computational region. The height should grow with streamwise distance to accommodate the growth in thickness of the boundary layer. The programs here are based on the untransformed equations to keep the codes simple, so the user must pick a computational region high enough at the initial station to accommodate the estimated growth of the boundary layer by the downstream end of the computational region. If not, the calculation must be stopped in the middle and then restarted with a higher region. The last matter concerns the choice of the step size(s). The integral methods only need a choice of dx, while the differential methods need dx and dy. This is normally accomplished by picking the number of points across the height and length of the region. For simplicity, most of the codes included here have fixed step sized in the x and y directions. For a laminar flow, about 20-25 points across the boundary layer are sufficient. Thus, knowing the thickness of the boundary layer at the initial station permits a reasonable choice for dy. With the boundary layer approximation, we can take dx >> dy, and a value of dx » d i/2 is usually adequate. Of course, for explicit methods, the step size must satisfy the stability criterion. For a turbulent flow a stretched grid in the y direction is much more efficient , but that was implemented in only one of the codes in this appendix in the interests of simplicity. Therefore, a large number of grid points is necessary, about 1000 points across the boundary layer.
This code implements a numerical solution method for Incompressible, Laminar Boundary Layer flows with an Implicit technique. See Schetz (1993) for details of the development. A default case corresponding to the following example problem is shown.
Example Laminar Implicit Numerical Method Problem: Consider 2D laminar flow of a fluid with a kinematic viscosityv = 2.0x10-4 m2/s at U¥ = 10.0 m/s over a surface that is a flat plate from the leading edge to x = 1.0 m. At that station, a ramp begins that produces an inviscid velocity distribution Ue(x) = 10.5 - x/2, m/s. This is an adverse pressure gradient, since Ue is decreasing so that p increases. Calculate the boundary layer development over this surface up to x = 2.0 m. Does the flow separate?
This is the same flow problem solved with the Thwaites-Walz integral method using the code WALZ and the explicit numerical method using code ILBLE. Now, we can apply the implicit numerical method in code ILBLI to this problem for comparison in terms of accuracy of the predictions and computational effort required.
Many, but not all, parts of the input are the same as for the explicit calculation. Since the first part of the problem is flow over a flat plate, the Blasius solution can be used to obtain "initial" conditions at x = 1.0. Thus, the numerical calculation will begin at XI = 1.0 and go to XF = 2.0. We must again provide input data for the kinematic viscosity as CNU = 0.0002 and the freestream velocity as UINF = 10.0. The Blasius solution at x = 1.0 gives the boundary layer thickness, DELT = 0.0224, the displacement thickness, DELTS = 0.00771, the momentum thickness, THETA = 0.00297 and the skin friction coefficient, CF = 0.00297. Choose 21 points across the initial boundary layer (DY = 0.00112) and add about 80 points above that to give MMAX = 100. Since the implicit method is unconditionally stable, no stability criterion need be followed, and we can select the streamwise step size DX and the number of points NMAX based only on accuracy considerations. Select NMAX = 41 to give DX = 0.025, which is about the size of the initial boundary layer thickness.
Initial profiles for U and V are required. For simplicity, adopt the cubic velocity profile for U and V = 0 as adequate approximations to the exact Blasius solution which would require tabular input for the profiles.
Finally, the inviscid velocity distribution is required. Since this is a linear edge velocity variation, velocities at only two points, at XI and XF need be specified. Again, UE is normalized with UINF in the code.
Press START and watch the skin friction, integral quantities and the u(x,y) velocity profile develop. Tabular values of the output can be accessed in the windows indicated.