Chinese Journal of Aeronautics, (2013),26(5): 1135-1146

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

cja@buaa.edu.cn www.sciencedirect.com

JOURNAL OF

AERONAUTICS

A parameterized geometry design method for inward turning inlet compatible waverider

Tian Chao, Li Ni *, Gong Guanghong, Su Zeya

Department of Automation Science and Electrical Engineering, Beihang University, Beijing 100191, China

Received 11 September 2012; revised 28 January 2013; accepted 13 March 2013 Available online 2 August 2013

KEYWORDS

Hypersonic inlets; Hypersonic vehicles; Inward turning inlet; Optimization; Parameterization; Waverider

Abstract Intensive studies have been carried out on generations of waverider geometry and hypersonic inlet geometry. However, integration efforts of waverider and related air-intake system are restricted majorly around the X43A-like or conical flow field induced configuration, which adopts mainly the two-dimensional air-breathing technology and limits the judicious visions of developing new aerodynamic profiles for hypersonic designers. A novel design approach for integrating the inward turning inlet with the traditional parameterized waverider is proposed. The proposed method is an alternative means to produce a compatible configuration by linking the off-the-shelf results on both traditional waverider techniques and inward turning inlet techniques. A series of geometry generations and optimization solutions is proposed to enhance the lift-to-drag ratio. A quantitative but efficient aerodynamic performance evaluation approach (the hypersonic flow panel method) with lower computational cost is employed to play the role of objective function for optimization purpose. The produced geometry compatibility with a computational fluid dynamics (CFD) solver is also verified for detailed flow field investigation. Optimization results and other numerical validations are obtained for the feasibility demonstration of the proposed method.

© 2013 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

1. Introduction

Intensive efforts have been devoted to the development of air-breathing hypersonic technology, where generation techniques of aerodynamic feasible geometries are accompanying the whole history. Due to the strong multidisciplinary cou-

pling nature and the high-velocity operating condition, the integration of an air-breathing system and a vehicle configuration is suggested to be a key step to achieve a favorable synthetic performance.

Three types of air-breathing solutions are conventionally adopted to configure hypersonic cruising vehicles, namely, the outward inlet (OI), the two-dimensional inlet (2DI), and the inward turning inlet (ITI). The ITI possesses the advantages of high pressure ratio, high mass flow capturing capacity, and robust off-design operation performances. Hypersonic vehicle configuration solutions are mainly categorized into wing-body type and lifting-body type, where the waverider configuration is fairly an important lifting-body option which is easy to be generated and managed to invoke intensive studies. However, the ITI system and

* Corresponding author. Tel.: +86 10 82339016. E-mail addresses: solarberiden@gmail.com (C. Tian), Lini@buaa. edu.cn (N. Li).

Peer review under responsibility of Editorial Committee of CJA.

1000-9361 © 2013 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. http://dx.doi.org/10.1016/j.cja.2013.07.003

waverider integration issues are rarely documented in the academic community probably due to complicated geometric properties and relatively expensive performance evaluation processes. Research of waverider and related air-breathing system integration and parameterized geometries is mainly restricted to 2DI-compatible geometries (X43A-like) or axi-symmetric conical flow field induced waveriders rather than free form surfaces based configurations. This situation limits the visions of developing new aerodynamic configurations for hypersonic vehicle designers. Nevertheless, a US patent1 shows that the NASA FALCON project claims to adopt the ITI system for the hypersonic air-breathing system while detailed integration approaches are unable to be clarified in public domains. Besides, it is known that the NASA FALCON project was later cancelled due to technical reasons and re-launched as the HTV-3X Blackswift project scheduled until the fiscal year 2023 from the official report of the US Department of Defense,2 which further points out the important trend of developing ITI-waverider integration techniques. The free form surface featured geometry design method proposed in this paper is easy to be implemented by off-the-shelf techniques and compromises three steps: (1) geometry parameterization; (2) performance evaluation (objective function); (3) performance optimization.

2. Overview of previous works

2.1. Geometry parameterization

The geometry parameterization step mainly focuses on the detailed means to integrate waverider design and ITI design.

On the aspect of waverider design, Sobieczky et al. proposed a hypersonic waverider design method from given shock waves3 and extended the method to adapt more complex parameterized supersonic cruising vehicles.4'5 An osculating flow field method and an osculating axisymmetric flow method of waverider generation led to more practical waverider configurations.6-8 These two methods could be extended with considering viscous effects9'10 and validated by experimental verifications.11

As for the topic of the ITI design, the approaches of devising three-dimensional hypersonic inlets with rectangular-to-elliptic section transition were detailed and summarized by Smart et al.12-16 The section transition inlet design techniques with intensive numerical simulated investigations were conducted by Gollan et al.17-19

Off-the-shelf works of designing the waverider and the ITI provide the prerequisites of the ITI-waverider compatible integration procedures. Therefore, it is urged to propose a novel method linking the two techniques. Fig. 1 shows the overall parameterization idea of the paper. The parameterized design method in this paper is extended by the inspiration of the HTV-3X configuration (see Fig. 1(a)). The optimization baseline is shown in Fig. 1(b). As the major surface is depicted by the non-uniform rational B-spline (NURBS) surface, the corresponding control points of the NURBS patches are shown in Fig. 1(c) and (d). The orthogonal sections in Fig. 1(e) illustrates the interior surface configuration.

Furthermore, a series of suggestions are made for obtaining complete hypersonic geometry generations since

(a) Original layout of NASA HTV-3X1

(b) Proposed layout in this paper

(c) The NURBS control points of the top surface

(d) The NURBS control points of the bottom surface

(e) The orthogonal sectioning illustration

Fig. 1 Overview of the hypersonic cruising vehicle geometry produced by the method proposed in this paper.1

few references could be listed to forge a complete CAE-ready-to-use aerodynamic configuration.

2.2. Performance evaluation and optimization

The aerodynamic performance is evaluated after the geometry parameterization step, which is intended to concentrate onto an objective function establishment for aerodynamic geometry optimization purpose. Considering the frequent calling of a performance evaluator, the classic hypersonic flow adapted panel method solver20'21 is utilized in the optimizer algorithm.

The waverider geometry optimization using differential evolution algorithm22-25 is aimed to enhance cruising efficiency gcrs. Cruising efficiency could be rephrased as the ratio of the specific cruising rang ds achieved to the specific fuel mass dM consumed according to the Breguet range equation26

= A = L 1 m

gcrs dM - ' n ' " ^

where dM is also the decremental variation of the cruiser's mass equaling the fuel mass consumption, Vc the cruising speed, L/D the lift-to-drag ratio, and c the specific fuel consumption. For a given propulsion system with a specific

airframe configuration, VcL/D needs to be maximized for enhancing the cruising performance. The scramjet powered waverider geometry optimization is practically recommended to complete the geometric integration procedures because the scramjet engine is difficult to efficiently heat the post-shock high-temperature gas flow to provide adequate thrust to cancel the drag.27 Cruising efficiency is therefore substantial on design specifications.

3. Detailed parameterization procedures

3.1. Overview of the generation procedures

The parameterization sequential procedures correspondingly numbered in Fig. 2 are listed as below:

(1) Generate the ITI geometry according to the on-design cruising Mach number of the vehicle.

(2) Extend from the ITI geometry to forge the complete dual-mode scramjet engine.

(3) Parameterize the leading edge curves and initialize the preset angle of attack (a0).

(4) Design the pseudo inlet capture curve and generate the leading edge surface which reflects the key procedures of the ITI-waverider integration method.

(5) Extend smoothly the downstream bottom surface with the NURBS surface following the leading edge surface.

(6) Complete the waverider top surface.

(7) Generate the downstream fairing system of gas exhausting system and vertical tails.

(8) Patch the geometric gaps for watertight features with constrained Delaunay triangulation methods if further computational fluid analysis is required.

All the curves and surfaces are generated by using the NURBS techniques. The above eight steps can render complete geometry data for both a computational fluid dynamics (CFD) solver and a hypersonic flow field adapted panel method solver.20,21

It needs to mention that all the dimensions in the current paper are not assigned with any unit systems. Use simply the 1-unit length for future engineering scaling treatment.

3.2. A brief review of previous related works

The ITI system is believed to be a promising air-breathing solution for hypersonic cruising technology, and the classic

ITI system is generated from a standard characterized flow field, which is normally obtained by three genres of equation solutions: the internal Taylor-MacColl equations,28 the method of characteristics equations,29,30 and the Navier-Stokes equations. After the standard characterized flow field is acquired, off-centered stream tracing and rectangular-to-ellipse transition12-19,31-33 techniques are offered to form the ITI geometry.

For a dual-mode scramjet engine, the configuration of an appropriately sized isolator34-36 and a diverged combination of a combustor and a nozzle by numerical simulation is believed to be necessary to enhance the stability and performance of the propulsion system.37 The integration of the ITI, the isolator, and the combustor and nozzle system is numerically tested to be able to generate positively installed thrust by duplicating the simulation configuration in Ref.37, which offers the complete engine system dimensions (see Fig. 3) for considering the overall hypersonic vehicle geometry generation. The scramjet engine component surfaces are symbolized by SSJ

(SITI, ^SOL SCONZ).

The traditional conical shock wave flow field induced wave-rider geometry generation approach stands on the assumption that the crosswise flow versus the mainstream flow is negligible. The streamwise flow traces settling on the osculating planes are sufficient to represent the three-dimensional flow fields. The conical flow induced waverider design method is based on a pure conical flow field by solving the Taylor-Mac-Coll equation (see (1) in Fig. 4), while the osculating conical flow induces waverider design method is an extended version of the pure conical one (see (2) in Fig. 4). The flexibility of the osculating plane philosophy to represent a variety of wave-rider configurations depends on the determination of the flow capture tube (FCT) and the inlet capture curve (ICC). Each osculating plane orientation is traditionally defined by the continuous and monotonous curvature normals along the ICC. The ICC is intentionally required to be designed on the same crosswise cutting plane for being installed on a scramjet and is initially devised for depicting the inlet mass flow capturing interface. Besides, from the practical point of view, the curvatures along the ICC can hardly be guaranteed to be adequately continuous and smooth for generating a stably induced bottom surface. These requirements for the ICC restrict the osculating cone waverider generation method.

3.3. The waverider leading edge

Since the FCT curve is essential to the conventional osculating conical flow field induced waverider generation, it is necessary to parameterize the waverider leading edge similar to the FCT,

Fig. 2 Schematic illustration of sequential procedures for a complete ITI-waverider integration geometry.

x(Spanwise^treamwise) O

Fig. 3 Overall dual-mode scramjet engine system dimensions.

Fig. 4 A brief review of conventional conical and osculating conical flow field induced waverider design principles.

which governs the general upwind profile of the produced vehicle. It is intended to reassemble part of the leading edge curve and part of the upwind edges of the ITI to produce the pseudo FCT (pFCT) (the dotted curve in Fig. 5). It should be noted to emphasize that the preset a0 significantly impacts the overall lift-to-drag ratio; hence an optimization process of a0 is believed to be necessary.

Given overall scramjet engine dimensions and a waverider leading edge curve, the key procedure of wrapping the ITI geometry with the exterior waverider geometry is resolved into a problem of finding an appropriate and less expensive approach to depict the near-corner transition surface between the ITI and waverider geometries using the osculating plane techniques (see Fig. 5).

3.4. Pseudo inlet capture curve

It is the essential philosophy of the osculating conical flow field method to assume that the crosswise flow is negligible and osculating plane slices of the streamwise flow is stackable. It is substantial to position and simplify the transitioning osculating planes' distribution along the pFCT to adapt the complicated near-corner flow field (see Fig. 5). As the essential counterpart of the pFCT, the ICC is renamed the pseudo ICC (pICC) to induce the leading edge surfaces (see Fig. 6) while the osculating plane is renamed the pseudo osculating plane (pOP) since the relation of the plane being tangential to the conventional ICC is dismissed in this paper. The leading edge surface induction is segmented into four types: normal

type, transition type, reverse type, and horizontal extrusion type (see Fig. 6). In Fig. 6, note that (1) some control points (such as the 8-9-10, 12-13, 3-4-5) are intentionally close to one another on purpose of NURBS curve generation; (2) the pFCT and the pICC are separated into counterparts of the Sections A, B, C, while Section D has no pFCT and pICC. The pFCT location could refer to the dotted curve in Fig. 5.

3.5. The horizontal extrusion surface

Domain D is a horizontal splaying extrusion of the ITI top lip curve and the curves on the opposite sides provide two segments of the pFCT. The splaying angles k1 and k2 are significant to eliminate the singular turning point near the corner of the flow field (see Fig. 6). The horizontal extrusion surface is symbolized by SHES.

3.6. The normal type induction surface

The pFCT and the pICC are practically represented by the same number of piecewise linearly discretized points for drawing corresponding pOP's. The normal type induction surface is generated by solving the conventional Taylor-MacColl equation. The line connected by the corresponding points AF (on the pFCT) and AI (on the pICC) denotes the line on the Mach cone surface (see Fig. 7), and the x-axis passing through AF denotes the axial direction of the Mach cone. This process simplifies the osculating plane method and facilitates the initialization of the Taylor-MacColl equation solver to generate a single curve on the induced bottom surface.

The coordinates of the pFCT and the pICC are sufficient to invoke the simplified pOP method. The Taylor-MacColl equation is rewritten as

(1 - v2 - vf)(2vr + vrcoth + v'0)- v'r(vrvr + W)= 0 (2)

where h is the angle variable between shock wave angle hs and hc, which are the shock cone angle and the solid cone surface angle (see in Fig. 7); c is the ratio of specific heat, and the derivatives are taken of h; vh is the angular velocity component and vr is the radical velocity component with the equation holding:

Vh = ^ = dh

Fig. 5 Waverider leading edge and pFCT parameterization and the preset angle of attack. Note that the separated control points LE0-LE10 define the NURBS leading edge curve and the dotted curve is the pFCT.

Fig. 6 pICC curve and the 21 control points (the pFCT and the pICC stretching up the bottom leading edge surface with three types of induction methods; Surface D being a simple horizontal splaying extrusion of the inlet top lip).

Fig. 7 Schematic illustration for the generalized pOP method.

In addition, the non-dimensional velocity V (with respect to the stagnation sonic speed) components are related to the local aft-shock Mach number Ma2 (see Fig. 7) by

v2 = v2 + vh

(y - 1)Ma\

Define

"X1 ' Vr (в) 'Vr'

X2 Х(в). .Ve.

Substitute Eq. (6) into Eq. (2) to yield a set of ordinary differential equations

x' = /(в, x) =

x2x1 — 1-1 (1 — Xl — x2)(2x1 + x2 cot в)

2 (1 Xj X2) X2

For the initial values of x and x2, the integration of Eq. (7) is started with the shock wave angle hs

hs = arc tan

|(ypICC — ypFCT)+ i(zpICC — zpFCT)|

XpICC — XpFCT

where XpKc > XpFcTypFcr = 0 and ZpFCT = 0 due to the locality of the origin of the local coordinates (see Fig. 7):

hs arc tan

lypICC + i • zpicc| xpICC _ xpFCT

kp = j(ypICc + i • zpicc)

where k means the argument function of the complex variable and kp is the argument of complex variable ypICC + i • zpICC. Therefore, the stream deflection angle behind the shock wave satisfies

ô = arc tan

2 cot hs

Ma1(sin2 hs — 1) Ma2(y + cos(2hs)) + 2

Normal components of the Mach number on the shock wave surface hold the following equations:

Man1 = Ma! sin hs (12)

Man2 =

Man,(y — 1)+2

2yMah(y — 1) — (y — 1)

where Man1 and Man2 are the normal components of the Ma! and Ma2 through the shock wave. According to Fig. 7 and Eq. (5),

V2 = 1

(y — 1)Ma^/ sin(hs — ô)

where V2 means the corresponding velocity to the Mach number Ma2, which is the Mach number after the local oblique shock. Hence the initial values of Eq. (7) are

r(0) _

x(0) x1 ' V2 cos(hs - d)

x(0) V sin(hs - d)

The integration of Eq. (7) is initialized by Eq. (15) until x2 = vh = 0 to get the velocity field (vr, vh) given any position in the local coordinates of the pOP (x, q) (see Fig. 7). The local velocity (vr, vh) and the local position (x, q) can be related by

/ : (x, q) # (vr, Vh) (16)

The normal type stream tracing integration begins locally with the point C:

- (xpFCT; l^pFCT + i ■ zpFCTj) - (xpFCT, 0)

until xpICC by solving the streamline equation

d(x, q)

- = /(x, q) = (vr, Vh)

Hence, On the pOP, the induced shock wave stream tracing curve can be piecewise linearly symbolized by a point sequence CjJ}n| forming a stream tracing curve:

C[0 ,n| -

xpFCT ■■■ xi xpICC

0 qi , ■■■, qn

The subscript i (i = 0,1,...,n) means the i-th point on the stream tracing curve on the local pOP coordinate (x, q). Convert the local pOP coordinate (x, q) to the three-dimensional coordinate (x, y, z) considering Eq. (10):

xpFCT xi xpICC

0 , ■■■ , q, cos Kp , ■■■, qn cos Kp

0 _ q, sin Kp _ _ qn sin Kp _

C|0 ,n] -

Combining all the sweeping curves Cj0 n]

of Domain A in

Fig. 6 produces the surface SN(CN n]). In this paper the Taylor-MacColl equation does integral by the lower order stiff differential equation solver ode23tb in MATLAB.

A robustness problem is addressed for the normal type induction surface generation: when the demanding integral shock wave angle hs drops below the critical angle of approximately hs c « 10°, due to stiffness and precision issues, varieties of solvers (including the multistep Runge-Kutta solvers, the multistep Runge-Kutta-Fehlberg solvers, the LSODE solvers from Lawrence Livermore National Laboratory, and the MATLAB ODE solvers) render reasonless results: surface wrinkling, integral warnings or divergence errors. The test results are listed in Table 1 when on-design flight condition is assigned with Mflra = 6, y =1.4.

This phenomenon can cause serious the surface generation failure when the integral demanding shock angle hs (hs < hs c)

Table 1 Failure shock wave angle hs c of solvers.

Solver name Approximate failure 0s c (°)

Multistep Runge-Kutta-Felberg 613.0

Multistep Runge-Kutta 613.1

MATLAB ode23tb 68.3

MATLAB ode45/ode113/15s 69.5

MATLAB ode23/ode23s 69.1

MATLAB ode23t 68.7

LLNL LSODE 612.2

Average failure 0s c 610.2

is assigned and can greatly influence the robustness of the induced surface generation procedure. An applicable spline interpolation remedy is prescribed for bridging the surface generation to enhance the robustness of the surface generation method:

(1) Collect a series of shock angles hs i(Ma1, c) and the corresponding cone angles hci(Ma1, c) and enforce the bridging function passing through the origin, while the bridging interpolation function is defined as

fbrg: [As,i(Mai,c) , 0|#[0c,;(Mai ,c) , 0]hB,i > hs,c - 10° (21) Given hs < 10°, the resulting cone angle is

hc -f(hs) , hs < 10° (22)

(2) Calculate the ending points qn in Eq. (19) with

(xpICC — xpFCT)

|qpICC qpFCTj

tan[hc(1 - r) + hsr]r -—0.25,hs < 10°

(3) Regulate the interpolation function fbrg in Eq. (21) if the on-design smooth condition is changed.

The bridging interpolation function is specified for a prescribed design condition (as Max and c are given) (see Fig. 8). Shifting to another design condition requires the rebuilding of the bridging interpolation function. Fifty times of tests about robustness are carried out with the qpICC descending to zero (as the demanding shock wave angle hs descending to zero as well) while the integrator of Eq. (7) failed, i.e., the ending point qn generation on the pOP is interrupted at the angle of hs c. On the design condition Mam = 6.0 and c = 1.4, the interpolation function could smoothly bridge the surface generation when the prescribed shock wave point vanishes. While changing any parameters of Mam and c can cause the surface wrinkling at the angle of hs,c since the fbrg has to be updated to the adaptation of new design condition.

Fig. 8 Mam and c impact the smoothness of bridging in scenario hs < hs c. Note that Ma^ = 6.0 and c =1.4 represent the on-design smooth condition, and other parameter configurations cause surface wrinkling with the pICC ending point decreasing.

Fig. 9 Normal type induction surface and the reverse type induction surface in scenarios of reflection weakening ITI and Busemann ITI.

3.7. The reverse type induction surface

The reverse type induction surface is located below the ITI system's lower cowl surface (see Fig. 6), where the wave-riding effect (caused by the upwind osculating conical flow field stream-tracing method) is not required while the drag reduction and downstream surface transition concerns come out to be important issues. Given a pICC, specifically near the ITI cowl tip section surface generation, the surface interference of the newly generated surface with the original ITI surface and the spillage drag reduction are two major factors to be considered. Using the reflection weakening ITI concept38 instead of the Busemann ITI design (see Fig. 9(b)), in order to alleviate the shock impingement and the reflection loss, enhance the pressure recovery ratio, and expand the surface interference margin of local surface generation, the reflection weakening ITI is intentionally designed to be downward curved (see Fig. 9(a)) with the streamwise trends after the impinging curved shock.

In Fig. 9, given a specific outer shock line, the original induced surface possesses a smaller upwind area than its mirrored counterpart. However, provided a specific upwind area, the original induction surface tends to interfere with the ITI interior wall especially near the cowl tip (see Fig. 9a). Despite in the Busemann ITI scenario, the original induction surface leads to a fairly thin upwind cowl lip, which may cause serious thermo protection issues. Furthermore, the original induction surface is designed to cope with the wave-riding features, the importance of which is replaced by the spillage drag reduction and surface interference considerations near the cowl tip section. Mirroring the normal induction surface about the Mach line into a reverse type is treated as an advantageous choice (see Fig. 9): both keeping the smooth continuity and utilizing the off-the-shelf algorithms. The reverse type induction surface is symbolized by all the curve set CR,[0 n] to produce the swept surface SR(CR [0 n]).

3.8. The transition type induction surface

Fig. 10 Schematic illustration of the transition type induction surfaces. Noting that T1 and S1 are the local curvilinear coordinate frame origins; the coordinate varies from T1/S1 to T2/S2; S0-S1 and T0-T1 are normal type induction surfaces; S2-T2 is a reverse type induction surface.

transitioned smoothly from the normal type surfaces (Domain A in Fig. 6) to the reverse type surfaces (Domain C in Fig. 6). The transition algorithm is illustrated in Fig. 10.

Taking T1-T2 local curvilinear coordinates as an example, during the transition procedures from T1 to T2, both the normal and reverse type induction surfaces are calculated along the positive direction of the curvilinear coordinates. The transition rate is controlled by an exponential factor n and generates the swept surface ST(Cj0 n), in which

Wi _ y^N

C [0,n] _ C [0,n]

-[0 , n]'

where the operation is imposed on the (x, y, z) coordinates of the sweeping curves CN0

"[0 ,n]> C[0 ,n]

T , and CR n], and

in which n = 4 is recommended. t (or s) is the curvilinear coordinate along the curve T1-T2 (S1-S2).

The bottom upwind surface SBU parameterization, consisting of SHES, SN, SR, and ST, is generally the most complicated and time-consuming part of the whole parameterization process. The resulted surface is illustrated in Fig. 11.

3.9. Downstream flexible surface generation

Downstream flexible surface generation is accomplished by carefully initializing the NURBS surface control points to form control curves in order to avoid the interferences between the scramjet geometry and the waverider outer skin surface patches. The initial NURBS control curves' configurations depend mainly on the designers' manual regulation. However, the downstream control curves are normally inherited from the corresponding upstream control curves for convenience. In order to keep the tangential relation with the leading edge surface, the downstream bottom surface starts with the linear extension of the rear edge of the leading edge surface (see Fig. 12). The top surface generation is similar to that of the

As illustrated in Fig. 6, the transition type induction surface is located in the Domain B, where the surfaces are required to be

Fig. 11 Resulted leading edge surface. Note that the dotted curve illustrates the pFCT curve.

Fig. 12 Bottom/top surface generation with the NURBS control curve.

Fig. 13 Overview of parameterization result.

bottom surface. The bottom downstream surface is named as SBD and the top surface as STP.

3.10. Parameterization summary

After the exhaust flow fairing geometries (named SFF) and the vertical tails (named SVT) are assembled, the parameterization process is seen in Fig. 13.

The parameterization result comprises 227 control variables and the major constraints for changing the variables are requested by the surface interference concerns. All the 227 control variables are listed in Table 2.

The geometry (see Fig. 13) is treated as a parameterized initial baseline values Wt, i = 1,2,.. .,227. The constrained range vector ei, i = 1,2,.. .,227 adding to each baseline value of W;-, i.e., (Wj + e) leads to geometry variations. An objective func-

tion is established for performance evaluation and optimization purpose. The resulted geometric surface data are symbolized by

G : (W, + e,) # S (26)

where i = 1,2,.. .,227, and S is a set of surfaces by:

S — fSSJ ; SBU ; SBD; STP ; SFF ,SVTg (27)

4. Performance evaluation

4.1. Overview of performance evaluation methods

Aerodynamic performance is a substantial design issue for hypersonic vehicles. Different levels of evaluation methods are appropriate for variety of design requirements. From the concerns of efficiency/precision trades-off, the conceptual design and optimization period leans towards computationally less expensive approaches which are feasible to effectively search a parameter space even with lower precision, and the panel method20'21 is one of such that is used by both industrial and academic practitioners. On the other hand, CFD solvers based on the finite volume method are intensively used to calculate detailed flow field parameters of interest. In this paper, it is necessary to demonstrate that the proposed geometry parameterization method is feasible to be evaluated by both the panel method and the finite volume based CFD method.

4.2. Parameterized geometry compatibility

The usual panel method requires simplified panel elements' representation of the waverider geometry, while the generated NURBS based geometry is fundamentally structured surfaces and naturally suitable for the panel element based aerodynamic performance evaluation programs. The panel elements' mesh and pressure coefficient Cp distribution results are illustrated in Fig. 14. It is notable that the panel method could give a rather fast performance evaluation by consuming approximately 0.3 s for each operating condition.

The preprocessing for the CFD solvers demands that the input geometry have to be watertight and un-sharp surfaces patched. Therefore, several blunting processes are necessarily implemented by the NURBS technique on the grid generation

Table 2 The 227 control variables summary.

i Control variables initial baseline position Wt Constraint range et Schematic figures

1-3 xLE0, УLE0, yLE1 [-10,10] Fig. 5

4 «0 [-1,2] Fig. 5

5,9,13,17 1 + |LE4 - LE3|/|LE3 - LE2| [-0.05,0.30] Fig. 5

1 + |LE6 - LE5|/|LE5 - LE4| [-0.01,0.30]

1 + |LE8 - LE7|/|LE7 - LE6| [-0.01,0.30]

1 + |LEm — LE9|/|LE9 — LE8| [-0.01,0.30]

6-8 xLE5, УLE5, zLE5 [-10,10] Fig. 5

10-12 xLE7, yLE7, zLE7 [-10,10] Fig. 5

14-16 xLE9, УLE9, zLE9 [-10,10] Fig. 5

18 xpICC0 [-5,5] Fig. 6

19-57 xpICC1,ypICC1v ■ ■,xpICC19,ypICC19, xpICC20 [-5,5] Fig. 6

58-77 zpICC0v . . zpICC20 [-5,5] Fig. 6

78-133 Bottom surface free control curve 1, 2, 3 [-5,5] Fig. 12

134-227 Top surface free control curve 1, 2, 3, 4, 5 [-5,5] Fig. 12

(a) Panel elements mesh of baseline geometry

(b) Pressure coefficient distribution of baseline geometry

Fig. 14 Initial baseline geometry adapted analysis by the panel element method.

(a) Unstructured grid generation by ANSYS-ICEM-CFD based on the watertight parameterized geometry

(b) CFD pressure contour result by ANSYS-FLUENT

Fig. 16 Compatibility demonstration of the parameterization geometry with grid generation software and a CFD solver.

concerns as well (see Fig. 15). The surfaces are exported with the stereo-lithography (STL) format.

For elaborate analysis of the parameterization result, the grid generator ANSYS-ICEM and the finite volume method based CFD solver ANSYS-FLUENT are utilized to validate the computability of the STL-formatted waverider geometry output. By configuring the upstream Mach number 6.8 and the hypersonic cruising altitude 30000 meters under fuel-off mode with parallel implicit time integrating density based solver and the k-e standard RNG turbulence model, one of the essential interests of using the CFD solver is to investigate the ITI system operating performance and the overall aerodynamic performance of the waverider. The rough computational result (see Fig. 16) intentionally demonstrates the watertight feature, the unstructured grid generation feasibility, and the computability of a complete CFD case based on the parameterized geometry and implies that the mesh generator program (ANSYS-ICEM-CFD) and the CFD solver (ANSYS-FLUENT) are compatible with the geometry.

4.3. Panel method based objective function

During the initial design phase, practitioners confront the issues of searching candidate geometry configurations. General qualitative and less precise quantitative performance data rather than massively detailed flow field parameters are of interest. Furthermore, a computationally inexpensive performance evaluator is provisionally the only means to close an optimization loop with an optimizer within a reasonable period of time since the optimization procedures intensively schedule the evaluator while heavy calculation loaded evaluators are normally unaffordable under present technological conditions.

For the hypersonic flow calculation issues in this paper, the panel method20,21 based hypersonic waverider performance evaluator (see Fig. 14) is employed to solve the lift-to-drag ratio as the objective function value. The objective function K is defined by

K: S # Cl/CD (28)

where S is provided by Eq. (26), while CL and CD are respectively the coefficients of lift and drag.

Fig. 15 The geometry preprocessing effects of blunting sharp edges using the NURBS technique. (1): The leading edge connecting the top and bottom surface smooth bridging; (2) and (3): joint spot smoothness of the bottom surface, the ITI surface, and the inlet top lip extrusion surface; (4),(8)-(11): the ITI surface and the bottom surface blunting spots; (5)-(7): smooth transition-ing between the leading and trailing edges of the vertical tails.

Waverider design requirements (Cruising Mach number, altitude,payload, ■ ,etc)

'V.+6 mete

v Geometry parameterization inputs

Waverider generation method

(¿VC*)

Performance evaluator(cost function)

Fig. 17 Parameterized geometry optimization control flow.

Table 3 The pseudo code of the DE geometry optimizer.

• Initialize the baseline geometry control values W, i = 1,2,.. .,227

• Set population number Np;

• Randomly initialize all Np agents ef, i = 1,2,.. .,227, j = 1,2,..Np in the population in the search space;

• Do until (maximum evolution generations produced)

• For each ef, j =1, Np

• Randomly pick three agents e(a),e(b), e(c) from the population, a,b,c 2 {1,2,...,Np}. e(a),e(b), e(c) are distinct from each other and all unequal to ef;

• Randomly pick R 2 {1,2,.. .,227};

• Compute the y, i = 1,2,.. .,227 potentially as the new optimal position for ef ;

• for each ef, i = 1,227

• Pick rt ~ U(0,1) uniformly from the open range (0,1);

• If (i = = R) or (ri < CR) then

(a) , !•/ (b) (c)\

• yi e\' + F( e)' — e) ');

• else

• y, := e,j);

• end if

• end for loop on i

• If (T(Wt + y.) > T(Wt+ e(j)) according to Eq. (29) then

• ef := yi;

• else

• continue loop on j

• end if

• end for loop on j

• Record the most improved sub-optimal agent e(j as the best found candidate solution; end do loop

Return the best solution e j) ever found.

5. Geometry optimization

5.1. Overview of geometry optimization

Geometry optimization is recommended for a complete scram-jet powered hypersonic waverider geometry generation since the air-breathing scramjet engine technically fails to efficiently heat the upstream gas after the hypersonic shock wave for generating adequate thrust. Geometry optimization is ultimately targeted to improve the overall cruising efficiency, which is impacted greatly by the lift-to-drag ratio.

The geometry optimization control flow is implemented as a computational loop (see Fig. 17). The initial waverider design requirements invoke the primary geometry parameterization inputs and activate the waverider geometry generation calculation. The calculation conveys the geometry information to the performance evaluator (i.e., the objective function K in Eq. (28)) which prepares the input for the optimizer. The optimizer closes the iteration loop to search for the optimal values from the baseline values for the parameterized waverider geometry within the constrained ranges.

5.2. The selection and configuration of optimizer

The objective function in this paper is naturally nonlinear and possesses 227 input variables (see Table 2). The differential evolution (DE) algorithm is selected to address such an optimization problem with nonlinearity and large input dimensionality. The DE optimization method is feasible to search large spaces of candidate solutions, but fails to guarantee that an optimal solution is ever found.

The basic searching mechanism of the DE algorithm functions by driving a population of candidate solutions

(called agents) to move in the search-space obeying simple mathematical formulae (variants of search schemes) to reference the existing agents from the population. The most improved agent is accepted as the sub-optimal solution in the current generation, and other agents are simply dropped and varied for next evolution generation. The objective function is practically the combination of Eqs. (26) and (28) leading to

T = G ° K: (Wi + ei) # CL/CD ( 29)

where W;-, i = 1,2,.. .,227, are the initial baseline positions of the control variables, and eh i = 1,2,.. .,227, are the baseline offsets to vary the resulting geometries. The ei offsets are randomly selected in the constrained range in Table 2. The pseudo code algorithm is listed in Table 3.

Note that F 2 [0,2] is called the differential weight factor and CR 2 [0,1] is called the crossover probability, and that the variants of the DE optimization schemes are consecutively developed to improve the optimization performance. In this paper, only the DE/rand/1 (see Table 1) variation scheme is utilized by the DE optimizer to verify the method's capability.

Table 4 Differential evolution (DE) optimizer settings.

Optimizer settings Option

Objective function dimension Bounded constraints Population number Np Evolution generations Differential weight factor F Crossover probability CR Differential strategy variants 227 ei constraints in Table 2 227 882 0.6497 0.9455 DE/rand/1

Fig. 18 Optimized waverider geometry.

100 200 300 400 500 600 700 800 Iterations

Fig. 19 CL/CD iteration history of the DE algorithm.

(a) Panel elements mesh of optimized geometry

(b) Pressure coefficient distribution of optimized geometry

Fig. 20 Panel element mesh and pressure coefficient distribution corresponding to the optimized waverider geometry.

The parameters Np, F, and CR as well as different variation schemes influence greatly the optimization performance. The optimizer settings of the present study are listed in Table 4. More recommended settings are summarized by Price et al.39

5.3. Optimization results and summary

The parameterized geometry design problem is implemented on a computer with CPU: Intel i7-2600/3.40 GHz and RAM: 16 GB. The iteration convergence history is illustrated in Fig. 18.

The optimized waverider geometry, panel element mesh, and pressure coefficient distribution contour map are shown in Fig. 18. The optimization result in Fig. 19 indicates that the DE optimizer is feasible to globally search the optimal solution that the lift-to-drag ratio has been elevated from less than 5.2 to 6.35, which implicates the enhancement of cruising performance.

A comparison between the panel method analysis results in Figs. 14 and 20 shows that the pressure coefficient distribution has been smoothly rearranged by the DE optimizer and demonstrates the necessity of the optimization in the parameterization geometry design procedures.

6. Conclusions

(1) The conventional osculating conical flow field induced waverider generation method is simplified and generalized to enhance the geometry flexibility and the optimization procedures which retain the effectiveness of waverider geometry generation.

(2) The leading surface generation is a combination of the induction surfaces with four types: normal type, reverse type, transition type, and horizontal extrusion type, which are the most time-consuming and complicated steps during the geometry parameterization phase.

(3) The applicability of the parameterized geometry has been tested with a commercial grid generator and a CFD solver for detailed flow field information. On the other hand, the less precise panel method has also been demonstrated to be capable to evaluate the pressure coefficient distribution of the waverider.

(4) The lift-to-drag ratio is selected as the objective function for the purpose of optimizing the parameterized wave-rider. Note that the optimization is recommended as a necessary phase for future scramjet powered waverider configuration design.

(5) The optimizer employs the differential evolution algorithm and accomplishes optimization results. The lift-to-drag ratio of the parameterized waverider in the present study has been elevated from less than 5.2 to 6.35, implying enhancement of cruising efficiency.

Acknowledgment

They would also like to thank the anonymous reviewers for their critical and constructive review of the manuscript. This study was supported by the National Natural Science Foundation of China (Grant No. 61004089).

References

1. Elvin JD, Lockheed-Martin Corporation. Integrated inward turning inlets and nozzles for hypersonic air vehicles. United States Patent US 7,866,599 B2. 2011, Jan 11.

2. Koch CMD. AF science, technology, and engineering overview. Headquarters US Air Force. <www.dtic.mil/ndia/2011SET/ Koch.pdf >; 2011.

3. Sobieczky H, Dougherty FC, Jones K. Hypersonic waverider design from given shock waves. In: Proceedings of the first international waverider symposium; 1990 Oct 17-19; 1990.

4. Sobieczky H, Choudhry SI, Eggers T. Parameterized supersonic transport configurations. In: Proceedings of the seventh European aerospace conference; 1994 Oct 25-27. Toulouse, France; 1994.

5. Sobieczky H, Stroeve JC. Generic supersonic and hypersonic, configurations. AIAA-91-3301; 1991.

6. Sobieczky H, Zores B, Wang Z, Qian YJ. High speed flow design using osculating axisymmetric flows. In: Proceedings of the third pacific international conference on aerospace science and technology; 1997 September 1-5. Xi'an, China; 1997.

7. Sobiecsky H, Zores B, Wang Z, Qian YJ. High speed flow design using the theory of osculating cones and axisymmetric flows. Chinese J Aeronaut 1999;12(1):1-8.

8. Qian YJ, Sobiecsky H, Eggers T. Waverider design with parametric flow quality control by inverse method of characteristics. In: Proceedings of international symposium on inverse problems in engineering mechanics; Nagano, Japan; 2000.

9. Geng YB, Liu H, Yao WX, Wang FM. Viscous optimized design of waverider derived from cone flow. Acta Aeronaut Astronaut Sin 2006;27(1):23-8 in Chinese.

10. Morgan MH, Khaikine VA. Waverider geometry generation using the osculating flowfield method with viscous effects. AIAA-2012-0951; 2012.

11. Miller RW, Argrow BM, Center KB, Brauckmann GJ, Rhode MN. Experimental verification of the osculating cones method for two waverider forebodies at Mach 4 and 6. AIAA-98-0682; 1998.

12. Smart MK. Design of three-dimensional hypersonic inlets with rectangular-to-elliptical shape transition. J Propul Power 1999;15(3):408-16.

13. Smart MK. Experimental testing of a hypersonic inlet with rectangular-to-elliptical shape transition. J Propul Power 2001;17(2):276-83.

14. Smart MK, White JA. Computational investigation of the performance and back-pressure limits of a hypersonic inlet. AIAA-2002-0508; 2002.

15. Smart MK, Trexler CA. Mach 4 performance of hypersonic inlet with rectangular-to-elliptical shape transition. J Propul Power 2004;20(2):288-93.

16. Smart MK, Ruf EG. Free-jet testing of a REST scramjet at offdesign conditions. AIAA-2006-2955; 2006.

17. Ferlemann PG, Gollan RJ. Parametric geometry, structured grid generation, and initial design study for REST-Class hypersonic inlets. In: 25th Propulsion systems hazards joint subcommittee meeting; La Jolla (CA), United States; 2009.

18. Gollan RJ, Smart MK. Design of modular, shape-transitioning inlets for a conical hypersonic vehicle. AIAA-2010-940; 2010.

19. Gollan RJ, Ferlemann PG. Investigation of REST-class hypersonic inlet designs. AIAA-2011-2254; 2011.

20. Gentry AE, Smyth DN, Oliver WR. The Mark IV supersonic-hypersonic arbitrary-body program, vol. I: User's manual. AFFDL-TR-73-159; 1973.

21. Gentry AE, Smyth DN, Oliver WR. The Mark IV supersonic-hypersonic arbitrary-body program, vol. II: Program formulation. AFFDL-TR-73-159; 1973.

22. Engelbrecht AP. Computational intelligence: an introduction. 2nd ed. San Francisco: John Wiley & Sons Ltd; 2007.

23. Karaboga D, Okdem S. A simple and global optimization algorithm for engineering problems: differential evolution. Turkey J Electric Eng 2004;12(1):53-60.

24. Storn R, Price K. Differential evolution-a .simple and efficient adaptive scheme for global optimization over continuous spaces. TR-95-012; 1995.

25. Storn R, Price K. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;ll(4):341-59.

26. Howe D. Aircraft conceptual design synthesis. London: Professional Engineering Publishing; 2000.

27. Segal C. The scramjet engine: processes and characteristics. New York: Cambridge University Press; 2009.

28. Van Wie D, Molder S. Applications of Busemann inlet designs for flight at hypersonic speeds. AIAA-92-1210; 1992.

29. Zucrow MJ, Hoffman JD. Gas dynamics, vol. I. New York: John Wiley & Sons, Inc.; 1976.

30. Zucrow MJ, Hoffman JD. Gas dynamics, vol. II: Multidimentional flows. New York: John Wiley & Sons, Inc.; 1977.

31. Billig FS, Baurle RA, Tam C, Wornom SF. Design and analysis of streamline traced hypersonic inlets. AIAA-99-4974; 1999.

32. Zhao Z, Song W. Effect of truncation on the performance of Busemann inlet. Mod Appl Sci 2009;3(2):168-71.

33. Sun B, Zhang KY, Wang CP, Jin ZG, Li N. Inviscid CFD analysis of hypersonic Busemann inlet. J Propul Power 2005;26(3):242-7.

34. Hutzel JR. Scramjet isolator modeling and control [dissertation]. Wright-Patterson Air Force Base (OH): Air Force Institute of Technology; 2011.

35. Smart MK. Scramjet isolators. RTO-EN-AVT-185; 2010.

36. Sullins G, McLafferty G. Experimental results of shock trains in rectangular ducts. AIAA-92-5103; 1992.

37. Milligan RT. Dual mode scramjet: a computational investigation on combustor design and operation [dissertation]. Dayton (OH):

Wright State University; 2009.

38. Zhang K, Nan X. Numerical and experimental investigation on hypersonic inward turning inlets with basic flowfield using arc tangent curve law of pressure rise. In: 17th AIAA international space planes and hypersonic systems and technologies conference; San Francisco, California; 2011.

39. Price KV, Storn RM, Lampinen JA. Differential evolution: a practical approach to global optimization. Berlin: Springer; 2005.

Tian Chao received his B.S. and M.S. degrees in Aircraft Propulsion and Power Engineering from Nanjing University of Aeronautics and Astronautics in 2008 and 2010, respectively, and has been a Ph.D. candidate in the School of Automation Science and Electrical Engineering at Beihang University since 2010. His main research interests are aerospacecraft geometry parameterization, high-efficiency aircraft performance evaluation, multidisciplinary design and optimization, computational fluid dynamics, etc.

Li Ni is an associate professor at Beihang University. She received both her B.S. and Ph.D. degrees in Control Science and Engineering from Beihang University. Her current research interests include modeling and simulation, decision support system and intelligent system and distributed computing technology. She is director & deputy secretary-general of Chinese Association for System Simulation (CASS) and secretary-general of Virtual Reality and Simulation Branch, Chinese Society for Stereology.