Sample Cover AIAA 2003-0774 Design-Oriented Aerodynamic Analysis for Supersonic Laminar Flow Wings I. Kroo, P. Sturdza Stanford University Stanford, CA 41st Aerospace Sciences Meeting & Exhibit 6 - 9 January 2003 Reno, Nevada For permission to copy or to republish, contact the American Institute of Aeronautics and Astronautics, 1801 Alexander Bell Drive, Suite 500, Reston, VA, 20191-4344. AIAA-2003-0774 Design-Oriented Aerodynamic Analysis for Supersonic Laminar Flow Wings Ilan Kroo∗ Peter Sturdza† Department of Aeronautics and Astronautics Stanford University, Stanford, CA Abstract Introduction The computation of boundary layer properties and laminar-to-turbulent transition location is a very complex problem, generally not undertaken in the context of aircraft design. Yet if an aircraft is to exploit the advantages of laminar flow and proper trade-offs between inviscid drag, structural weight, and skin friction are to be made, this is just what must be done. This paper summarizes a designoriented method for the aerodynamic analysis of supersonic wings including approximate means for estimating transition and total drag. The boundary layer analyses employed here are meant to be inexpensive but sufficiently accurate to provide some guidance for advanced design studies or to be incorporated in a multidisciplinary optimization. The boundary layer analysis consists of a quasi3D finite difference method, based on inputs from a 3D Euler analysis. A parametric fit of integral boundary layer properties to the amplification rate predicted by linear stability analysis, leads to a simplified method for predicting TS transition and the onset of instability of stationary crossflow modes. This simplifed boundary layer analysis and inviscid CFD were combined to evaluate the drag of a simple wing planform, defined by approximately ten geometric parameters. These parameters were used as design variables in a numerical optimization procedure to minimize total wing drag, illustrating the trade-offs between skin friction and inviscid drag. Three dimensional Computational Fluid Dynamics (CFD) has become a useful and popular tool for the analysis and design of aircraft. However, for multidisciplinary optimization (MDO) at the conceptual design stage, for cases in which laminar to turbulent transition play an important role in configuration drag, it is still far too expensive. Theory and limited experimental evidence suggest the possibility of maintaining extensive natural laminar flow over wings with moderate sweep at supersonic speeds. The wing planform, twist, and airfoil profiles must carefully designed, taking into account the effects on both aerodynamics and other disciplines. The tradeoffs between different disciplines are a key element to the viability of the concept. For example, wings that support the maximum extent of natural laminar flow may have little sweep and thus require very small thickness to reduce wave drag. This adversely affects the structural design, fuel volume, and subsonic performance. A detailed and complete multidisciplinary design study is necessary to determine whether these tradeoffs result in a more efficient configuration than a conventional design. (For a discussion of the MDO problem, see Manning and Kroo.1 ) Since 3D CFD is incapable of the quick execution times necessary for this type of design problem, a less expensive, yet reasonably accurate method of boundary layer computation within an optimization framework is sought. Such a method must provide a calculation of skin friction, including sufficient modeling to enable boundary layer transition prediction, while demanding little computational resources, since the analysis will likely be called thousands of times in the course of a nonlinear search procedure. In the following sections, we describe this approach to the design of supersonic wings with lam- ∗ Professor, † Graduate Fellow AIAA Student c Copyright 2003 by the authors. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission. 1 American Institute of Aeronautics and Astronautics 2 inar flow, including boundary layer calculation and approximate transition prediction. Comparisons with more refined calculations and application to a simple wing design problem illustrate the use of the method. Direct Optimization Approach Previous laminar design studies have often relied on the specification of a target pressure distribution, followed by nonlinear inverse design to determine a geometry compatible with laminar flow. The desired Cp may be defined based on the experience of aerodynamics experts and designers who understand some of the basic trade-offs between inviscid and viscous drag as well as the need for maintaining certain structural properties. However, recent experiences with the design of low boom, small supersonic aircraft suggest that this approach quickly becomes intractable when complex trades between supersonic and subsonic performance, boom signature shape, and wing structural weight are required. Since wing skin friction may account for less than 20% of the drag of a business jet design at Mach 1.6, the implications of the target Cp on properties other than laminar flow extent may be paramount. Furthermore, if only the outer wing panel is amenable to extensive laminar flow and/or it is only practical to exploit laminar flow on the upper (or lower) surface, the importance of laminar flow is further reduced and it is not clear that the best design should even try to obtain low skin friction at the expense of other characteristics. A quantifiable evaluation of this trade-off is possible using a direct optimization approach that involves varying the wing geometric properties and computing the effect on overall aircraft performance and constraints. This multidisciplinary optimization approach addresses the difficulties in identifying suitable target pressures and can provide important insight into basic configuration choices. It has not, however, been used for advanced design of laminar aircraft for several reasons. The computation of boundary layer properties and transition is extremely computation intensive and often not reliable. Results may also be rather noisy, causing potential difficulties with conventional search methods. The goal of the present research was to develop methods, which although necessarily approximate, would provide more direct guidance in the design of supersonic vehicles employing laminar flow for improved performance. Sweep/Taper Boundary Layer Analysis The approach to estimating boundary layer properties is based on an approximate sweep/taper method described by the authors in.2 Some of the discussion here is excerpted from that paper, while a more detailed description of the computations are provided in Sturdza’s thesis.3 The sweep/taper method was originally suggested by Bradshaw as an extension to the ideas used for calculating boundary layers on infinite yawed wings. The sweep/taper simplification allows the computations to be two dimensional even though a three dimensional flow field is solved. It applies to semi-infinite straight tapered, swept wings with straight isobars parallel to the spanwise generator lines. Bradshaw et al.,4, 5 apply this simplification to their turbulent boundary layer code and report that it runs at approximately two thirds the speed of their compressible 2D code. It does not compute laminar boundary layers due to the special numerical procedure of the Bradshaw-Ferriss-Atwell finite difference method. Kaups and Cebeci6 developed a laminar finite difference boundary layer program based on similar assumptions, however, it can only calculate the boundary layer up to the chordwise station where the crossflow changes sign (unacceptable for the current application). In the mideighties, Ashill and Smith7 developed an extension to Smith’s general 3D integral method that takes advantage of the Bradshaw sweep/taper theory. This approach, however, introduces additional assumptions that also make it unsuitable for this research. The main difficulty is, again, an inability to handle configurations where the crossflow changes sign. The present method is capable of laminar and turbulent calculations with a slightly improved version of Bradshaw’s formulation of the sweep/taper theory. The basic sweep/taper theory can be developed as follows. In the case of the infinite (untapered) yawed wing for incompressible laminar flow, it is well known that the boundary layer equations become two dimensional in the plane normal to the spanwise coordinate, with a third, decoupled equation for the spanwise velocity (Prandtl’s indepen- 3 ∂u ∂v + =0 ∂x ∂y ∂ ∂ ∂p ∂ ∂u x-mom: ρ u +v u=− + µ ∂x ∂y ∂x ∂y ∂y ∂ ∂ ∂ ∂w z-mom: ρ u +v w= µ , ∂x ∂y ∂y ∂y U 8 dence principle): cont: where (u, v, w) and (x, y, z) are the cartesian velocity components and coordinates in the streamwise, normal, and spanwise directions, respectively. In the compressible or turbulent case, the equations become coupled, but since the derivatives in the spanwise direction are zero, the equations are still computationaly two dimensional (there are only two independent variables). After adding taper to the infinite yawed wing, as discovered by Bradshaw et al.,4 the boundary layer equations can be transformed in order to remain computationaly two dimensional (with the additional assumption that the surface isobars are parallel to the spanwise generators of the wing). The 3D compressible boundary layer equations are written in cylindrical coordinates along an arc from the leading edge to the trailing edge that is perpendicular to both edges and to the isobars in between (Figure 1): cont: x-mom: z-mom: temp: ∂ρu ∂ρv ∂ρw ρw + + − =0 ∂x ∂y ∂z r o ∂ ∂ ∂ ρuw ρ u +v +w u− ∂x ∂y ∂z ro ∂p ∂ ∂u =− + µ ∂x ∂y ∂y ∂ ∂ ∂ ρu2 ρ u +v +w w+ ∂x ∂y ∂z ro ∂p ∂ ∂w =− + µ ∂z ∂y ∂y ∂ ∂ ∂ ρcp u +v +w T ∂x ∂y ∂z ∂p ∂p ∂ ∂T =u +w + κ ∂x ∂z ∂y ∂y " 2 # 2 ∂u ∂w +µ + . ∂y ∂y The origin of the coordinate system is where the wing would come to a point if it were semi-infinite in span. In this coordinate system, there is no pressure gradient in the spanwise, z, direction, and other derivatives in the z direction can be written in terms of derivatives in the direction normal to the wing z r0 Λ x θ Figure 1: Cylindrical coordinate system used in the sweep/taper boundary layer calculation (view of wing planform from above). surface which results in a set of computationally 2D equations. This comes from noting that along a generator line (constant x), the boundary layer thickness is closely proportional to r0 in turbulent flow,4 √ and r0 in laminar flow; essentially a conical flow assumption where the boundary layer can be viewed as self similar in the spanwise direction. Therefore the spanwise derivative, ∂f /∂z, of any flow quantity, f , (other than pressure) can be re-written in terms of a vertical derivative of that same flow quantity, (y/r0 )∂f /∂y for turbulent flow or (y/2r0 )∂f /∂y for laminar flow. These operations on the boundary layer equations result in numerically two dimensional equations that can be solved on a 2D x − y grid along an arc of constant radius from the origin of the coordinate system. The Bradshaw sweep/taper assumption limits, strictly speaking, the isobars to be parallel with the spanwise generators and thus constrains the lift distribution to a triangle. For the semi-infinite tapered wing, calculations along a single arc will provide the 3D flow solution on the entire wing. However, by running calculations at several spanwise stations and including some corrections to the basic theory described below, this approach can be used as a quasi3D calculation on a finite wing, taking into account the local effects of sweep and taper as long as the isobars are reasonably straight and parallel with the generators. Fortunately, it is not uncommon for the isobars to be somewhat parallel with the generator lines in both subsonic and supersonic wings of medium to high aspect ratio. The present code is an extension of the Bradshaw sweep/taper theory. Bradshaw’s method, as has been implemented previously, only uses freestream conditions and surface pressures along the computational arc as inputs. The external (inviscid) veloc- 4 ity components are calculated from those inputs using the sweep/taper assumptions. It has been found that the method can be improved by using more information from the inviscid solution. The entire inviscid flow field is available as an external boundary condition to the sweep/taper code. For example, the pressure gradient in the generator direction, ∂p ∂z , which is neglected in the original method, is calculated by the inviscid flow solver and can be applied directly in the boundary layer equations. The external velocity components are similarly available, but the difficulty is that they are not necessarily consistent with the sweep/taper equations. That deviation, however, can be used to develop corrections to the sweep/taper version of the mass, momentum and energy equations in the form of source terms. The present sweep/taper method implements these improvements by solving the following set of equations for laminar flow: ∂ρu ∂ρˆ v 3ρw + − = ρQ ∂x ∂y 2ro ∂ ∂ ρuw x-mom: ρ u + vˆ u− ∂x ∂y ro ∂ ∂u ∂p + µ + Fx =− ∂x ∂y ∂y ∂ ρu2 ∂ z-mom: ρ u + vˆ w+ ∂x ∂y ro ∂p ∂ ∂w + µ + Fz =− ∂z ∂y ∂y ∂ ∂ temp: ρcp u + vˆ T ∂x ∂y ∂p ∂ ∂T ∂p +w + κ =u ∂x ∂z ∂y ∂y " 2 # 2 ∂u ∂w +µ + + B, ∂y ∂y cont: where, vˆ = v + w y , 2ro and, ∂ρw y ∂ρw − ∂z 2ro ∂y ∂u y ∂u Fx = −ρw − ∂z 2ro ∂y ∂w y ∂w Fz = −ρw − ∂z 2ro ∂y ∂T y ∂T B = −ρcp w − . ∂z 2ro ∂y Q=− 1 ρ (1) (2) (3) (4) Mach alt. ΛLE ΛT E b/2 cr ct t/c (ft.) (deg.) (deg.) (ft.) (ft.) (ft.) (%) Agrawal 2.0 40,000 32.0 -15.4 18.94 22.5 5.45 5.0 NLF 1.5 45,000 16.7 -29.0 23.15 25.3 5.52 1.75 Flight Test 1.4 10,000 19.5 -27.5 3.0 4.16 1.57 2.5 Table 1: Geometry and flight conditions for the calculations appearing in this paper. The source terms, equations 1-4, and the spanwise pressure gradient, ∂p ∂z , are the corrections applied to the Bradshaw sweep/taper theory to arrive at the present formulation. The pressure gradient is independent of y, and is therefore known throughout the boundary layer, but the source terms are not. At the edge of the boundary layer, the source terms can be evaluated using the inviscid flow solution from an Euler or panel code. At the wing surface, w and y are zero, so the terms vanish there. In this implementation, they are modeled linearly in between.∗ These developments significantly improve sweep/taper results without any discernible computational cost over the original Bradshaw formulation of the theory. The sweep/taper equations are discretized and solved using Keller and Cebeci’s algorithm for the compressible boundary layer equations.8 In brief, the code solves the continuity, x- and z-momentum, and temperature equations in cylindrical coordinates along the circular arc drawn in Figure 1. The continuity equation is integrated out with stream functions and the remaining equations are transformed via a compressible Falkner-Skan procedure using the arc length and the tangential external velocity component in the similarity variables. The discretization is second order, and Newton’s method is used to march the solution along the arc from leading edge to trailing edge. Theoretically, the only limit on wing sweep and taper is that the marching always be downwind. Figure 2 presents the crossflow Reynolds number predicted using the sweep/taper method compared to detailed 3D Navier-Stokes computations (CFL3D)9 on the Agrawal wing planform (Table 1). Inviscid pressures and velocities used for all the boundary layer work in this paper were performed ∗ Actually, they are linear in the transformed Falkner-Skan vertical coordinate. 5 1 1 Present Sweep/Taper 3D Navier-Stokes Present Sweep/Taper 3D Navier-Stokes 0.8 0.6 0.6 y/δ Sweep/Taper without corrections Sweep/Taper with corrections 3D Navier-Stokes 1000 0.8 y/δ 1200 800 0.4 0.4 0.2 0.2 0 -0.025 Rcf 0 0 600 0.2 0.4 0.6 0.8 1 u/ue 400 0.6 y/δ y/δ 0.6 0.4 0.4 0.2 0.4 0.6 0.8 0 0.8 0.8 0 -0.005 Present Sweep/Taper 3D Navier-Stokes Present Sweep/Taper 3D Navier-Stokes 0 -0.015 -0.01 w/ue 1 1 200 -0.02 1 X/C 0.2 0.2 0 0 0.2 0.4 0.6 u/ue Figure 2: Crossflow Reynolds number: sweep/taper vs. 3D Navier-Stokes on the Agrawal wing at midsemispan. Laminar 3D Navier-Stokes Euler + Sweep/Taper Skin Friction 0.00043 0.00043 Inviscid 0.01923 0.01919 Turbulent 3D Navier-Stokes Euler + Sweep/Taper Skin Friction 0.00362 0.00356 Inviscid 0.01931 0.01919 Table 2: Wing drag coefficients: Sweep/taper vs. Navier-Stokes on the Agrawal wing for 100% laminar and 100% turbulent computations. with the FLO87 Euler code developed by Jameson.10 The present sweep/taper calculations are in quite good agreement with Navier-Stokes, showing the dramatic improvement over the original Bradshaw formulation. Not shown are Agrawal’s results11 which cannot go beyond the zero-sweep chord location because he used a modification of the Kaups and Cebeci code. The velocity profiles calculated near mid-chord and near the trailing edge with the present sweep/taper method also compare favorably with Navier-Stokes results (Figure 3). It has been found that subtleties in the flow are captured at least qualitatively by this method. For instance, variations at different spanwise stations (such as wing root and tip effects—departures from strict sweep/taper assumptions) are computed accurately. This is shown in Figure 4 on an example NLF planform (Table 1), verifying that the sweep/taper method can be used as a “quasi-3D” theory by per- 0.8 1 0 -0.006 0 0.006 0.012 w/ue Figure 3: Streamwise (left) and crossflow (right) velocity profiles at 56% chord (top) and 92% chord (bottom) calculated with sweep/taper theory (solid line) vs. 3D Navier-Stokes (symbols) on the Agrawal wing. forming calculations at various spanwise locations. The coefficient of friction in the freestream direction is compared between the sweep/taper method and 3D Navier-Stokes in Figure 5 at one spanwise station. Complete wing drag coefficients are shown in Table 2 for laminar and turbulent computations on the Agrawal wing. The agreement is very good especially in the laminar flow case. Since different turbulence models are involved (Spalart-Almuras in CFL3D and Cebeci-Smith in the sweep/taper code) the difference of less than one count in turbulent skin friction is very reasonable. The comparisons confirm that the sweep/taper approach can accurately calculate friction drag on thin supersonic wings and that there is no need for viscous-inviscid iteration at these flight conditions. Figure 6 depicts the temperature profile in the boundary layer calculated with the sweep/taper method and with Navier-Stokes. The two sets of Navier-Stokes data are for the same run, but at two convergence levels (at overall residual levels of 10−6 and 10−8 ). It takes an order of magnitude more iterations to achieve temperature convergence compared with the ususal CL or CD convergence criteria — two processor-weeks on a 195 MHz SGI Origin in this case. This behavior has consistently been observed in 2D Navier-Stokes calculations as well.12 6 400 350 Sweep/Taper Outboard Sweep/Taper Inboard 3D Navier−Stokes Outboard 3D Navier−Stokes Inboard 10-2 SWPTPRBL CFL3D 250 200 150 Cf Crossflow Reynolds Number 300 10-3 100 50 0 0 0.1 0.2 0.3 0.4 0.5 X/C 0.6 0.7 0.8 0.9 1 10-4 0 Figure 4: Crossflow Reynolds number: sweep/taper (lines) vs. 3D Navier-Stokes (symbols) at two spanwise stations on an NLF wing. Since temperature strongly influences the stability of Tollmien-Schlichting waves, full Navier-Stokes solutions must be highly converged to be useful in stability calculations. For this reason, the speed advantage of the present sweep/taper method over 3D Navier-Stokes is quite large. A week or two of runtime can be shortened to about 12 seconds on the same computer. Transition Prediction Transition to turbulence is one of the most complex aspects of fluid dynamics, however, reasonable semiempirical methods exist that can capture the general behavior of the transition front. One of the most common of such methods is linear stability theory and the en transition criteria. The en method is the basis of the current work. As additional 3D supersonic experimental data become available, improved techniques can be substituted in the general framework developed here. Transition to turbulence in a boundary layer is caused by growing disturbances in the laminar boundary layer. These flow disturbances are initiated mainly by surface roughness or undulations, mechanical and acoustic vibrations and by shock waves or turbulence in the free-stream flow. The disturbances start out as small, periodic perturbations to the mean flow, then grow or decay depending upon the characteristics of the boundary-layer velocity and temperature profiles. At first, when the disturbances are small, their behavior is similar to sinusoidal plane waves. As the unstable modes grow 0.2 0.4 0.6 0.8 1 X/C Figure 5: Skin friction coefficient: Sweep/taper vs. Navier-Stokes on the Agrawal wing. in amplitude, nonlinear interactions become dominant. Following the nonlinear growth, these laminar disturbances begin to spurt out intermittent spots of turbulence which spread and eventually merge together resulting in a fully turbulent boundary layer. A good discussion of the physics of transition can be found in White.13 These laminar instabilities are classified into several different types including Tollmien-Schlichting waves (TS), G¨ortler vortices, attachment line instabilities, and crossflow vortices (CF). The focus of the present work is on the development of approximate methods by which TS and CF transition can be estimated. Transition Analysis Overview A complete analysis of transition can really only be done with direct numerical simulation (DNS), however, as in the computation of turbulence, this method is limited to small Reynolds numbers due to the enormous computational requirements. Even so, knowledge of the factors that originally excite these instabilities is incomplete — the initial or boundary conditions for the DNS. In the case of an airplane wing, it is difficult to specify a complete picture of the conditions: free-stream turbulence, acoustic energy both in the flow and in the structure, microscopic details of the surface shape and surface finish and contaminants. Receptivity, the response of laminar modes to various disturbances, is an active area 7 1 0.8 Present Sweep/Taper 3D Navier-Stokes (Resid: 10-8) 3D Navier-Stokes (Resid: 10-6) y/δ 0.6 0.4 0.2 0 200 220 240 260 280 300 320 340 360 T (Kelvin) Figure 6: Temperature profile showing slow convergence of N-S calculation (the more converged case required an order of magnitude more iterations). of research. Although rules of thumb and empirical results from receptivity analysis and experiments can be very useful, detailed receptivity analysis is not reasonable at the conceptual and preliminary design stages. The solution of stability equations represents a step down from the brute force DNS approach. These equations are a form of the Navier-Stokes equations after subtracting out the steady mean flow terms, leaving just the description of the unsteady disturbances. The significant advantage is that with relatively few assumptions, these equations can be manipulated to allow very efficient numerical solutions. Since one of the approximations results in partial differential equations with parabolic character, the equations are known as the parabolized stability equations (PSE). The simplifying assumptions mainly involve a description of the perturbations to the mean flow as a periodic phenomena characterized by a slowly varying amplitude function and rapidly varying wave behavior. The basic idea is that the oscillatory portions are transformed by Fourier series, leaving just the amplitude function to solve. Since the amplitude, or shape function, varies slowly in the streamwise direction, much less computational resolution is required to capture an accurate solution compared with DNS solutions. Additionally, a careful choice for the decomposition of the solution into the wave and shape functions allows some terms to be neglected, resulting in parabolic equations that can be solved with very efficient space marching schemes. Further details can be found in papers by Malik, Chang and Li.14, 15 The accuracy of the solution is excellent because laminar instabilities conform very nicely to the simplifying assumptions in the PSE, especially the nonlinear form of the PSE.16 Although the PSE provide a practical method with today’s computers to calculate laminar instabilities, there are drawbacks. For transition caused by streamwise Tollmien-Schlichting waves it is widely accepted that PSE is not necessary since Linear Stability Theory (LST, described below) provides good transition prediction capabilities. For the other dominant mode in swept wings in atmospheric flight – stationary crossflow vortices – it is thought that LST is adequate if the transition amplitude is empirically found for similar geometries.17 Additionally, it is not simple to predict transition for the stationary crossflow mode using PSE because the crossflow vortices themselves do not cause transition directly, but rather saturate and allow a secondary instability to grow which then triggers transition. This means that yet another stability code is used to calculate the growth of this secondary instability and an nfactor correlation is then applied. Another drawback is that the nonlinear PSE require as inputs initial conditions describing the birth of the laminar instabilities which are not generally available. Frequently, when detailed experimental flow fields are not known, the PSE solutions are started from mode shapes derived using LST. LST, on the other hand, is very well suited for this work. It can be used as a self-contained, though approximate, transition prediction method without requiring extensive knowledge of initial conditions. The main drawback is that the numerical procedures are poorly-suited to optimization. This is because LST solvers need to be run interactively, with an expert user monitoring results, discarding non-physical solutions, modifying inputs, and helping the solver along in some cases. Otherwise, LST solvers are reasonably fast, can accommodate the high Reynolds number flows in this study and do not require special starting procedures. It is true that LST does not provide as detailed a description of the physical flow disturbances as PSE, however that is not always an advantage for transition prediction at the conceptual design stage. Greater resolution of the flow requires a deeper understanding of the physics involved in the process of receptivity and transition itself. The simpler theory can be more robust. 8 LST can be summarized as follows. The NavierStokes equations are the starting point. They are separated into mean flow terms and fluctuation terms which are then simplified by neglecting nonlinear effects. Then the mean flow equations are subtracted out, leaving linearized differential equations governing the small disturbances. The parallel flow assumption is applied, which is basically a statement that the boundary-layer growth is small over a wavelength of the disturbance. With the additional assumption of sinusoidal solutions (via complex exponentials), the equations simplify into an eigenvalue problem at each chordwise station. This means that the solution procedure involves marching from leading to trailing edge and solving for the local disturbance properties. The inputs are the disturbance frequency and wave angle, and the outputs are the amplification rate and mode shapes for those inputs. Mathematical detail can be found in a paper by L. M. Mack.18 The difficulty is in the numerical procedure that finds unstable modes and tracks them as they travel downstream. This quality makes linear stability theory (LST) codes necessarily interactive, with the user providing reasonable frequencies and checking for non-physical results, or modes that are lost from station to station. For this reason, even a very powerful computer will still not enable the use of LST codes in numerical design optimization. Linear stability theory, and for that matter PSE calculations as well, don’t predict transition, but rather just predict the growth of periodic laminar instabilities. The most popular transition criterion based on LST calculations originated in the mid 1950s. Smith and Gamberoni, and, independently, van Ingen compared linear stability calculations with experimental results of various boundary layers where the transition location was known.19, 20 They found that the transition location in many experiments coincided with the location where linear stability theory predicts the amplification of TollmienSchlichting waves to reach about 8100. Since the commonly used amplification n factor is the natural logarithm of the amplification ratio, this became known as the e9 method. This method works very well for two-dimensional boundary layers that experience transition caused by growth of Tollmien-Schlichting instabilities. It is also assumed that the free-stream turbulence and acoustic disturbances are not too great since a long, slow, linear growth of the instabilities is required for LST to be accurate. If the initiating disturbances are large, then nonlinear growth can dominate the transitional region and violate the assumptions of LST. Low turbulence tunnels and atmospheric flight conditions generally meet these requirements. The n factor can be adjusted to other than 9 to account for varying receptivity effects, such as free-stream turbulence or surface roughness. L. M. Mack suggested this approach by proposing a correlation for n as a function of free-stream turbulence. He argued that the free-stream turbulence influences the initial amplitude of the laminar instabilities.21 As a rough guide, transition in wind tunnels generally observed at n factors of 7 to 9 whereas in atmospheric flight, it is more common to see transition at n factors of 10 to 14. For this reason, this criterion is now more commonly known as the en method where n is found empirically for the conditions of interest. Present Transition Prediction Methodology The present development assumes that 3D transition to turbulence on swept wings can be split into streamwise and crossflow criteria that are calculated separately and subsequently recombined. Compressible linear stability theory is heavily used since it successfully predicts TS dominated transition, and, it is believed, stationary crossflow transition on swept wings. Since the goal is to provide a transition analysis that can be used for aircraft design optimization, linear stability results are predetermined and stored as parametric fits. These fits are evaluated during the optimization to calculate instability amplification. The en criterion is used to predict transition for both the streamwise TS calculation and the crossflow modes (but with different values of n). Streamwise Instability Calculation Figure 7: Boundary-layer instabilities for similarity solutions. The streamwise methodology is borrowed from the work of Drela,22 Gleyzes23 and Stock and Dagen- 9 hart.? Most of these previous works are only applicable to incompressible flows, with Drela extending the method to the transonic regime. 8 7 6 Present Fit 2500 Hz 2750 Hz 3000 Hz 3250 Hz 3500 Hz 4000 Hz 4500 Hz 5000 Hz 5500 Hz 6000 Hz 7000 Hz 8000 Hz 9000 Hz 10000 Hz 12000 Hz 15000 Hz 17500 Hz 20000 Hz 25000 Hz 30000 Hz 5 15 N M=1.5 t/c=1% 4 Amplification Ratio, n 3 2 M=1.5 t/c=2% 10 1 M=2.4 t/c=1% M=1.5 t/c=3% 0 0.5 1 1.5 2 2.5 x (m) 5 M=2.4 t/c=3% 0 0 0.2 0.4 0.6 Chordwise Station, x/c 0.8 1 Figure 8: Amplification ratios of TollmienSchlichting waves on biconvex airfoils: present method (solid lines) compared to COSAL. 9 8 7 6 Present Fit 2500 Hz 2750 Hz 3000 Hz 3250 Hz 3500 Hz 4000 Hz 4500 Hz 5000 Hz 5500 Hz 6000 Hz 7000 Hz 8000 Hz 9000 Hz 10000 Hz 12000 Hz 15000 Hz 17500 Hz 20000 Hz 25000 Hz 30000 Hz N 5 4 3 2 1 0 0.5 1 1.5 2 2.5 x (m) Figure 9: Amplification ratios of TollmienSchlichting waves a biconvex airfoil: present method compared to LASTRAC. The fundamental idea is that a table lookup or regression curve fit can be made by parameterizing the instability growth rate. For the incompressible case, similarity solutions are used to relate growth rate to the boundary-layer shape factor and momentum thickness Reynolds number. Figure 7 is an example of how the rate of increase of n, the logarithmic amplification ratio, can be approximated by a linear envelope over a range of unstable frequencies. Since the Tollmien-Schlichting waves don’t amplify below a certain Reynolds number, that quantity is also modeled as a function of shape factor. The end result is a calculation of n that corresponds to the most unstable frequency at that chordwise station, however the value of the frequency is not found. The present approach builds upon Drela’s para- Figure 10: Amplification ratios of TollmienSchlichting waves on a flight test airfoil: present method compared to LASTRAC. metric formulas for two-dimensional amplification rates as functions of boundary-layer shape factor, H, and momentum thickness Reynolds number, Reθ . He actually used the so-called kinematic versions, Hk and Reθk , claiming that these formulations are less sensitive to compressibility effects up to transonic free-stream velocities. The kinematic formulas are just the incompressible versions of the H and Reθ equations where the density is not inside the integrand. The argument is that Tollmien-Schlichting growth rates depend much more on velocity profiles shapes than the density-weighted compressible parameters. For the current work, Drela’s formulation is modified in two ways. First, it was originally designed around his integral boundary layer method, so many of the empirical formulas for determining Hk and Reθk can be dispensed with since those quantities are calculated directly. Second, and more importantly, his method does not account for the stabilizing effect of a supersonic free stream, so a new parameter is introduced. The ratio of wall temperature to the external, inviscid temperature is added as the third parameter. The current method uses terms with this new parameter to modify Drela’s original formulas. It is based on linear stability results on biconvex airfoils at Mach 1.5 to 2.4. The parametric fits of amplification rate and critical Reynolds number are expressed in the following manner: dn Tw (5) = f Reθk , Hk , dReθ Te Tw Reθ0 = f Hk , . (6) Te As the boundary-layer equations are solved by the usual downstream marching, Reθk , Hk and Tw /Te 10 are evaluated at each station. The TS amplification factor can then be determined by integrating Equation 5 from the point at which the flow exceeds the critical Reynolds number defined in Equation 6 to the current chordwise station: Reθ Z nts = on transition prediction. The next higher step in fidelity from the crossflow Reynolds number method is, as for the Tollmien-Schlichting case, parametric fits of linear stability results. crossflow Reynolds number = dn dReθ . dReθ crossflow shape factor = Reθ0 Figures 8 to 10 show n-factors calculated with the present method compared with linear stability results. Crossflow Instability Calculation Crossflow transition is more difficult to model than the streamwise TS modes. The physical phenomena is predominantly nonlinear and requires, for greatest accuracy, nonlinear PSE solutions and an additional secondary instability code. Presently, only a few people in the world can accomplish this for subsonic flows, much less on a supersonic aircraft. Consequently, the methods used during conceptual and preliminary design of aircraft have been very crude to date. Prior to the present work, the most widely used design-oriented crossflow transition criterion was based solely on the crossflow Reynolds number, defined as: ρe |wmax |δcf , Rcf = µe where δcf , the crossflow thickness, is the height above the surface where the crossflow velocity component is 1/10th of its maximum, wmax . The method referred to as the “crossflow Reynolds number criterion” in this work is due to Malik et al.24 who introduced a correction for compressibility and stated that transition is predicted when γ−1 2 Rcf = 200 1 + Me . 2 This formula is based on correlations to experiments and computations on cones and swept wings in supersonic flow. It unfortuately turns out that the crossflow Reynolds number approach is often too conservative for transition prediction on a typical natural laminar flow wing. The primary problem is an over sensitivity to wing root and tip disturbances and to wing sweep. Since low wing sweep has an adverse effect on inviscid drag and sonic boom, it is important to be as realistic as possible in determining sweep effects crossflow velocity ratio = ρew max δcf µe ymax δcf w max Ue δcf w max y max Figure 11: Crossflow profile and parameters used for simplified crossflow instability growth model. The parameters used for the stationary crossflow fit are summarized in figure 11. These parameters were suggested by Ray Dagenhart25 and used by him to successfully model growth rates calculated with an incompressible linear stability code. Crossflow n-factors are calculated by modeling the amplification rates and integrating them while simultaneously solving for the boundary layer flow: wmax , Hcf (7) αδcf = f Rcf , Ue Rcf0 = f (Me ) Zx ncf = αdx. (8) (9) x0 Hcf is the crossflow shape factor and α denotes the spatial amplification rate of the instabilities. The point at which the stationary crossflow modes become unstable, x0 , is determined with a crossflow Reynolds number criterion using Malik’s compressibility correction, γ−1 2 Rcf0 = 43 1 + Me . 2 The value of 43 works quite well from Mach 1.4 to 2.2 on a number of different airfoils at several sweep angles. The previous parametric models for incompressible Tollmien-Schlichting instabilities were developed using boundary-layer similarity solutions. Unfortunately, there are no similarity solutions possible for compressible 3D boundary layers. Instead, to help identify the influence of each parameter on 11 Figure 12: Crossflow velocity profiles used in “similarity sequence” boundary layers: baseline, adjusted to increase crossflow shape factor, and adjusted to decrease crossflow velocity ratio. Model of stationary crossflow amplification rate at Mach 1.8 4 3.5 Amplification Rate, −α i 3 2.5 2 1.5 1 0.5 0 −0.5 0 0.2 0.4 0.6 0.8 1 X/C Figure 13: Stationary crossflow amplification rate for “similarity sequences” showing LASTRAC data (symbols) and the parametric model. Model of stationary crossflow amplification ratios at Mach 1.8 7 6 5 N for Hcf N for basic N for cfvr model for Hcf model for basic model for cfvr N 4 3 2 1 0 0 0.2 0.4 0.6 0.8 1 X/C Figure 14: Crossflow n-factors calculated with LASTRAC (symbols) and the parametric fit. Figure 15: Crossflow n-factors calculated with LASTRAC (symbols) and the parametric fit on three airfoils from Mach 1.4 to 2.2. growth rates, “similarity sequences” are used. A boundary layer is constructed by repeating the same profile (streamwise and crossflow velocities and temperature) at each chordwise x station and √ just scaling the dimensions in proportion with x. Therefore, the boundary layer is made of a sequence of profiles with constant crossflow shape factor and velocity ratio. This way each parameter can be individually adjusted to capture its contribution to the magnitude of instability amplification. It is expected that such an approach is reasonable when linear stability codes are used on the similarity sequence boundary layers because LST eigenvalue solutions are local to each x station. If PSE were utilized, however, this type of boundary layer may result in erroneous solutions because the conservation of mass, momentum and energy is violated in the mean flow. An example is presented in Figures 12-14. Figure 12 shows the crossflow velocity profiles used for three similarity sequences: the baseline and modifications that independently vary the crossflow shape factor and crossflow velocity ratio. Figure 13 presents the resulting amplification rate of the stationary crossflow modes. The linear stability solver in LASTRAC26, 27 is used to calculate amplification rates, and a parametric fit corresponding to Equation 7 is devised to follow this data. In this case, an analytic funtion is found that approximately fits the LASTRAC data, but a table-lookup method could also be used. Finally, Figure 14 shows the resulting integrated n-factors. It turns out that the three parameters picked by Dagenhart are sufficient to model crossflow distur- 12 modeled with a composite amplification ratio, q 2 2 N ∗ = (nts /ntscrit ) + (ncf /ncfcrit ) . The present transition criteria is preliminary, but it is hoped that for the cases that apply to laminar supersonic aircraft, the current approach will suffice. Nevertheless, correlations with a full three dimensional en stability code28 should be performed to validate and improve the streamwise and crossflow calculations used here, as well as the method used to combine them into a 3D criteria. Figure 16: Comparison of stationary crossflow amplification between COSAL and the present parametric method on a tapered wing. bances in supersonic flows. Figure 15 shows surprising accuracy for this example parametric fit that is based on data at only one Mach number. This observation confirms the assertion by Dagenhart that crossflow instabilities are much less influenced by compressibility compared to the TS modes. Since LASTRAC cannot calculate flows on tapered wings, the present fit is developed with only infinite-swept wing flows. It can be seen in Figure 16 that the fit compares favorably to COSAL28 results on a tapered wing with 30 degrees of aft leading edge sweep and 15 degrees of forward sweep at the trailing edge (the small-scale test blade flown on an F-1529 ). 3D Transition Criteria The interference between TS and crossflow instabilities is still being debated. Dagenhart25 briefly discusses the issue. He cites Pfenninger’s previous work and says that as long as both TS and crossflow modes are not simultaneously strongly amplified, their interaction can be ignored. Indeed, when the disturbances are small in amplitude and accurately described by linear theory there would be no interference. But Dagenhart does state that, “relatively weak oblique TS waves can distort and stretch crossflow disturbance vortices to produce rapid, resonancelike amplification and transition.” The current approach to this problem is to maintain simplicity because it is not well understood. Two methods have been used. In one, the crossflow and TS modes are taken as completely separate, with their own critical n-factor (9 for TS and 5 for CF). In the second method, a slight interference is Example Design Application As an example application of the boundary layer and transition estimation methods, a simple wing geometry was optimized to identify the minimum total drag design and to illustrate the trade-off between inviscid and viscous drag components. The example also served to evaluate the feasibility of optimization tools that include viscous drag and transition modeling for more complex wing design work in the future. By employing a direct optimization approach, the trade-offs between pressure gradient effects on TS and CF transition and on inviscid wave drag and structural constraints may be incorporated in a rational way. The software used in this exercise consisted of a number of computer codes including an inviscid Euler solver, boundary layer computations, transition and skin friction calculations, a response surface generation routine, and a nonlinear search method, integrated into an automated system and executed on an inexpensive cluster of personal computers. The wing inviscid (Euler) analysis was completed using Jameson’s FLO107 with a refined grid for accurate drag and inputs to the boundary layer analysis. A parallel (MPI) version of this code was used to provide reasonable turnaround times. Next, boundary layer properties were computed based on the Euler results and the new sweep/taper analysis code. Properties of the boundary layer were computed along 12 arcs across planform. From this data, cross-flow and TS transition were estimated based on the previously described parametric fit of LASTRAC data using a combined transition estimate based on an interaction-free model. Finally the skin friction drag is estimated using results from the laminar analysis and a Cebeci-Smith model for the turbulent boundary layer which is added to the inviscid drag from the Euler code. 13 Figure 17: Section geometry variations based on simple parameterization. In this study a response surface / trust region optimization method was employed to avoid potential problems with non-smooth viscous results. This approach involves generating quadratic models of the objective function and constraints, then applying a constrained nonlinear optimization procedure to identify the optimal result based on the fit. A revised fit is then computed centered around the new estimated optimum point and the process repeated. Hundreds of cases were analyzed to obtain multiple quadratic fits in the course of the optimization, providing detailed boundary layer data in each case for 12 sections over span. The ten design variables included angle of attack and twist along with several section properties at the root and tip stations. A linear variation of properties between root and tip was assumed. Section geometry variations included changes in thicknessto-chord ratio, location of maximum thickness, a parameter that affected forward and aft curvature, and a parabolic camber amplitude. Figure 17 illustrates the variations in thickness distribution represented with these degrees of freedom. The objective was to minimize total drag while constraining the overall wing CL and maintaining a constant section moment of inertia at root and tip to maintain torsional stiffness. Figures 18- 19 include selected results for this example wing with a 40 deg leading edge sweep at Mach 1.6. In Figure 18, the optimized wing’s amplification factors for TS and cross-flow are depicted. This wing, which yields minimum total drag, achieves almost full chord laminar flow over much of the upper surface, while obtaining only 14% laminar flow on the lower surface in a compromise to maintain section torsional moment of inertia and maintain low inviscid drag. The optimizer achieves the large upper surface laminar extent by reducing the curvature on the forward portion of the sections, carefully trading TS Figure 18: TS and cross-flow amplification factors on upper surface of optimized wing. Figure 19: Typical section Cp distribution for optimized wing. and CF growth. The lower surface is very stable to TS, but would need some type of control for crossflow stability on the forward part of section. The Cp distribution, shown in Figure 19 shows the very linear upper surface gradient over the full span. This exercise was repeated to define the wing geometry with minimum inviscid drag. Not surprisingly, this calculation, skipping the boundary layer analysis, led to a wing with essentially biconvex sections and a slightly higher inviscid L/D than the wing optimized for minimum total drag. When the viscous analysis was performed on this optimal inviscid wing some laminar flow was still predicted, but the total (viscous) L/D was estimated to be 8% lower. 14 Conclusions A design-oriented aerodynamic analysis was developed to permit direct optimization of supersonic wings employing laminar flow. This method, based on a modified sweep/taper concept and fits to linear stability theory results, permits rapid computation of boundary layer properties including estimates of transition and total drag. Although approximate, it provides important guidance in the advanced design of laminar wings. The boundary layer analysis, coupled with an Euler solver and numerical optimization method, was used to design a wing with minimum total drag at Mach 1.6. Extensive natural laminar flow was possible on this wing even with relatively high leading edge sweep of 40 degrees. Future work will focus on improving transition correlations, more efficiently handling the very nonlinear optimization problem, and applying the method to a variety of design problems. Acknowledgments The authors would like to thank Jeffrey Viken of Innovative Aerodynamic Technologies, Inc. for providing some of the en results; Chau-Lyan Chang at NASA Langley, author of LASTRAC who provided help with his code and general advice on transition prediction; David Lednicer at AMI who generated Euler solutions of the F-15 test wing, and Peter Bradshaw who supplied sweep/taper programs, documentation and answered many questions. Discussions with Profs. William Saric and Helen Reed of Arizona State University, Prof. Eli Reshotko, and Drs. Peter Hartwich and Art Powell of Boeing also have been especially helpful. These studies were initially and continually inspired by our interaction with Dr. Richard Tracy and ASSET Research Corporation, interested in the application of laminar flow to efficient supersonic aircraft. Much of this work has been supported by Directed Technologies, Inc. and their critical role in this research is gratefully acknowledged. References [1] Manning, V. M., and Kroo, I. M., “Multidisciplinary Optimization of a Natural Laminar Flow Supersonic Aircraft,” AIAA Paper No. 993102, Proceedings of the 17th Applied Aerodynamics Conference, 28 Jun. 1999. [2] Sturdza, P., Manning, V. M., Kroo, I., Tracy, R., “Boundary Layer Calculations for Preliminary Design of Wings in Supersonic Flow,” AIAA Paper No. 99-3104, Proceedings of the 17th Applied Aerodynamics Conference, 28 Jun. 1999. [3] Sturdza, P., “Design Optimization of Supersonic Natural Laminar Flow Aircraft Including Transition Prediction,” Ph.D. Dissertation, Stanford University, 2003. [4] Bradshaw, P., Mizner, G. A., and Unsworth, K., “Calculation of Compressible Turbulent Boundary Layers on Straight-Tapered Swept Wings,” AIAA Journal, Vol. 14, Mar. 1976, pp. 399-400. [5] Bradshaw, P., Unsworth, K., “An Improved Fortran Program for the Bradshaw-FerrissAtwell Method of Calculating Turbulent Shear Layers,” Imperial College Aero Report 74-02, Feb. 1974 (Revised Nov. 1980). [6] Kaups, K, and Cebeci, T., “Compressible Laminar Boundary Layers with Suction on Swept and Tapered Wings,” Journal of Aircraft, Vol. 14, Jul. 1977, pp. 661-667. [7] Ashill, P. R., and Smith, P. D., “An Integral Method for Calculating the Effects on Turbulent Boundary-Layer Development of Sweep and Taper,” Aeronautical Journal, Vol. 89, Feb. 1985, pp. 43-54. [8] Cebeci, T., and Bradshaw, P., Physical and Computational Aspects of Convective Heat Transfer, Springer-Verlag, New York, 1984. [9] Rumsey, C., Sanetrik, M., Biedron, K., Melson, N., and Parlette, E., “Efficiency and Accuracy of Time-Accurate Turbulent Navier-Stokes Computations,” Computers and Fluids, Vol. 25, No. 2, 1996, pp. 217-236. [10] Jameson, A., and Alonso, J., “Automatic Aerodynamic Optimization on Distributed Memory Architectures,” AIAA paper 96-0409, 34th Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 1996. [11] Agrawal, Shreekant, and Powell, Arthur G., “Supersonic Boundary-Layer Stability Analysis of an Aircraft Wing,” Journal of Aircraft, Vol. 28, Nov. 1991, pp. 721-727. 15 [12] Goodsell, A. M., “Computational Analysis of Blunt, Thin Airfoil Sections at Supersonic and Subsonic Speeds,” Ph.D. Dissertation, Stanford University, May 2002. [24] Malik, M. R., Balakumar, P., and Chang, C.L., “Linear Stability of Hypersonic Boundary Layers,” Paper No. 189, 10th National AeroSpace Plane Symposium, 23-26 Apr. 1991. [13] White, Frank M., Viscous Fluid Flow, 2nd ed., McGraw-Hill, New York, 1991. [25] Dagenhart, J. R., “Amplified Crossflow Disturbances in the Laminar Boundary Layer on Swept Wings with Suction,” NASA TP 1902, Nov. 1981. [14] Chang, C.-L., Malik, M. R., “Non-Parallel Stability of Compressible Boundary Layers,” AIAA paper 93-2912, 24th Fluid Dynamics Conference, Orlando, FL, July 1993. [15] Malik, M. R., Li, F. and Chang, C.-L., “Crossflow Disturbances in Three-Dimensional Boundary Layers: Nonlinear Development, Wave Interaction and Secondary Instability,” J. Fluid Mech., Vol. 268, 1994, pp. 1-36. [16] Haynes, T. S., and Reed, H. L., “Simulation of Swept-Wing Vortices using Nonlinear Parabolized Stability Equations,” J. Fluid Mech., Vol. 405, 2000, pp. 325-349. [17] Saric, W., personal communication, Arizona State University, Tempe AZ, 2001. [18] Mack, L. M., “Boundary-Layer Linear Stability Theory,” AGARD Rept. 709, 1984, pp. 3-1 to 3-81. [19] van Ingen, J. L., “A Suggested Semi-Empirical Method for the Calculation of the Boundary Later Transition Region,” Tech. Rep. VTH-74, Inst. of Tech., Delft, 1956. [20] Smith, A. M. O. and Gamberoni, N., “Transition, Pressure Gradient, and Stability Theory,” Douglas Aircraft Report ES-26388, 1956, also Proc. Ninth Internat. Cong. Appl. Mech., Vol. 4, 1957, pp. 234-244. [21] Mack, L. M., “Transition Prediction and Linear Stability Theory,” AGARD Conf. Proceedings 224, Paris, 1977, pp. 1-1 to 1-22. [22] Drela, M., “Two-Dimensional Transonic Aerodynamic Design and Analysis Using the Euler Equations,” Ph.D. Dissertation, MIT GTL Rept. No. 187, Feb. 1986. [23] Gleyzes, G., Cousteix, J., and Bonnet, J. L., “A Calculation Method of Leading-Edge Separation Bubbles,” in Numerical and Physical Aspects of Aerodynamic Turbulent Flows II, edited by T. Cebeci, Springer-Verlag, New York 1984. [26] Chang, C.-L., “The Langley Stability and Transition Analysis Codes (LASTRAC) – Part I: LST, Linear & Nonlinear PSE for 2D, Axisymmetric and Infinite Swept Wing Boundary Layers,” NASA TM, in preparation, 2003. [27] Chang, C.-L., “The Langley Stability and Transition Analysis Codes (LASTRAC): LST, Linear & Nonlinear PSE for 2D, Axisymmetric and Infinite Swept Wing Boundary Layers,” AIAA paper 2003-0974, 41st Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 2003. [28] Malik, M. R., “COSAL—A Black Box Compressible Stability Analysis Code for Transition Prediction in Three-Dimensional Boundary Layers,” NASA CR-165925, May 1982. [29] Kroo, I. M., et al., “Natural Laminar Flow for Quiet and Efficient Supersonic Aircraft,” AIAA paper 2002-0146, 40th Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 2002.

© Copyright 2020