Assessment of methods for laminar wing boundary layer calculation

The present study examines the methods available for calculation of laminar boundary layers on wings in subsonic/transonic flow, with an eye toward including such calculations within an existing optimizer-riven design algorithm. Since the stability and transition of the laminar flow is to be assessed, both the solutions and their smoothness are of interest. Robustness, speed of execution, and ease of selecting input parameters are also important. The performance of a number of computer programs was compared using three test cases: a flat plate, an airfoil, and a swept wing. The programs showed only slight differences in the calculated values and smoothness of boundary layer characteristics, but large differences in robustness, execution speed, and ease of selecting input parameters. Natural or suction-stabilized laminar boundary layers can yield reduced viscous drag on wings for transport aircraft. The present study aimed to evaluate available computer programs for laminar boundary layer calculation for their possible use in optimizer-driven wing design. To assess stability and transition behavior, detailed calculation of the laminar wing boundary layer is required. Many computer programs have been developed in the last 20 years which are capable of solving some form of the equations which govern the development of a compressible laminar boundary layer. It is now also possible to solve the complete Navier-Stokes equations without the boundary layer approximation but this approach was * Research Assistant, Member AIAA 7 Associate Professor of Mechanical Engineering, Member AIAA $ Aerospace Engineer, Advanced Aerodynamic Concepts Branch, Member AIAA not investigated here, as the repeated solutions required for wing design would be too costly in computer time. A number of computer programs were acquired for the present study; those to be discussed here are listed in Table 1. We are interested in laminar boundary layers developing on wings in transonic (compressible) flow, but two programs which solve the 2-D, incompressible equations (one based on Thwaites' method, the other on the Cebeci-Keller approach1) are included for comparison purposes. The mathematically simplest set of equations which may be solved are the boundary layer forms of the 2-D conservation of mass, momentum, and energy equations for a compressible fluid. These equations approximately model flow over an unswept wing, originating at a stagnation point. A further approximation the independence principle may be invoked to use the 2-D equations to model the flow normal to a swept wing's leading edge. Programs I B L ~ , ~ , STAIWSTAN~~, and VGBLP~ solve the 2-D boundary layer equations. They treat the fluid properties and the energy equation somewhat differently, and employ different numerical solution approaches as well (see Table 1). The next step in mathematical complexity of the governing equations would be to model a swept wing as infinite, with no gradients in geometry or pressure parallel to the wing leading edge. Conicallysymmetric flow is a generalization of this, allowing for a tapered geometry with no gradients along the generating rays of the wing. Programs WING6.7 and BLSTA* solve the conical flow equations, while BL3D-9o9>l0 solves the infinite sweep equations as well as those for a fully 3-D boundary layer. Copyright O 1994 by the American Institute of Aeronautics and Astronautics, Inc. No copyright is asserted in the United States under Title 17, U. S. Code The U.S. Government has a royalty-fiee license to exercise all rights under the copyright claimed herein for Governmental purposes. All other rights are reserved by the copyright owner. 35 7 Solution of the full 3-D boundary layer equations represents the highest level of mathematical complexity for the boundary layer approach. In addition to requiring solution of a more complex set of equations, the 3-D model requires treatment of conditions at the wing root and tip. For example, the root plane may be approximated as a symmetry plane, or infinite wing conditions could be assumed. The tip region is often ignored i.e., the boundary layer equations are solved only up to the edge of the tip's zone of influence. Calculations of the flow over the wing proceed from the root flow and solution of the attachment line flow. The only program we acquired which is capable of solving the 3-D boundary layer equations is B L ~ D ~ O ~ ? ' ~ . However, the test cases presented here do not employ this program's 3-D solution capability. The programs compared here vary not only in respect to the equations solved, but also as to initial and boundary condition treatment and the numerical solution approach employed. Also, the selection of grid parameters influences the solutions obtained. The approach taken here has been to try to employ "comparable" density in conditions and solution grids, with perhaps somewhat more detail than might be selected for "production" use of the programs. A comparison of results for three test cases will be presented to summarize the performance of the programs. All test cases represent laminar, subsonic, compressible flow with an adiabatic boundary and no transpiration. The test cases are summarized in Table 2, which shows that the three cases represent successive increases in the level of complexity of the boundary layer flow. First, a flat plate calculation was selected to confirm basic trends with respect to the effects of compressibility on the 2-D boundary layer. Next, an airfoil test was selected to isolate the added effect of pressure gradient on each program. The third case is a swept wing, for which the boundary layer flow is three-dimensional. The first test case is simply a boundary layer developing on an adiabatic flat at a Mach number of 0.8 and a unit Reynolds number of 5 million per meter. Freestream air properties for this calculation are taken as those of the U.S. Standard Atmosphere at an altitude of 12.2 km. The second case is an unswept, lifting airfoil with a shape designation ~ S ~ ~ ~ ( 1 ) 0 2 1 3 l . This shape was designed to achieve natural laminar flow under cruise conditions with moderate load. We have calculated the pressure distribution for this airfoil using the computer program ~ ~ 0 6 s ~ ' ~ ; both the pressure distribution and the airfoil shape are shown in Figure 1. The FL06SD calculation employed 128 points around the airfoil. This mesh was deemed suEicient to provide a meshindependent solution by noting that negligible differences are obtained for a 256-point calculation. The test case only uses the upper surface results over the initial 60% chord. The third case involves a went w i ~ g designed for natural laminar flow at cruise Mach numbers near 0.8. The leading edge sweep of this wing is 25.3 degrees, and the wing geometry includes taper, twist, and a planform break. The pressure distribution for this wing is taken from a recent wind tunnel test; Figure 2 shows the pressure dstribution on the wing upper surface at the 40% semispan location. The experimental data has been interpolated to provide a more dense specification of freestream velocity or pressure coefficient for input to the boundary layer calculations, but not smoothed. Figure 3 shows all the calculated results for the flat plate test case, for a region extending 1 meter. The influence of compressibility is apparent: the displacement thickness is increased about 15% from the incompressible value. Momentum thickness remains largely unaltered from the incompressible values, so the boundary layer shape factor H increases about 15% also: from 2.6 for the incompressible flat plate, to about 3.0. The skin friction is, for all the programs, also very little changed from the incompressible values. One program, STAN7, requires a starting velocity profile to be provided to start the calculation, as it cannot start from a similarity solution. This was provided very near the leading edge using the Blasius solution. For the other two test cases, this characteristic makes STAN7 inconvenient to use; other programs are able to start from the stagnation point or attachment line. For this reason STAN7 results were not obtained for the other two cases. Somewhat unusual treatment was required to obtain the flat plate solution from WING, which solves the conical flow equations and must start from an attachment line. Very small sweep and taper angles (0.01 deg.) were required, and the boundary conditions actually specified a smoothly-accelerating flow near an attachment line, followed by a region of constant pressure. It was found that WING was very sensitive to the spacing and density of the specified boundary conditions. Many trials with varying spacing and density of the boundary conditions in the accelerated flow region were required to produce a solution for this case. Some oscillation of the solution is noted. Another note concerns the program IBL. This program requires that wall temperature be specified as a boundary condition; it is incapable of handling a specified (zero) heat flux boundary condition except in an approximate way, appropriate only for unity Prandtl number. The solution procedure implied by this characteristic would be repeated trials, varying the specified wall temperature to achieve the desired adiabatic solution. In practice, however, we found that sgeclfylng a wall temperature equal to that computed from Buseman's relation with a recovery factor equal to &' r provided very closely adiabatic solutions.