Applied Mathematics and Computation 174 (2006) 1539–1563
www.elsevier.com/locate/amc
Scramjet missile design using genetic algorithms R.J. Hartfield *, J.E. Burkhalter, R.M. Jenkins Aerospace Engineering, Auburn University, 211 Aerospace Engineering Building, Auburn, AL 36849, United States
Abstract The objective of this effort was to show that scramjet powered vehicle designs which are optimized for a given flight condition can be found using predictive modeling tools and a genetic algorithm (GA). This required the assimilation of basic modeling tools for the external aerodynamics, the engine internal ballistics for thrust and additive drag predictions, and for the mass and inertia properties. The GA routines for optimizing the design have been used in other design optimization efforts. The basic models have been assembled and specific designs have been developed. Results indicate that the flight envelope for a scramjet is very narrow and dimensions are critical if the engine is to develop thrust at a given design point. For the designs considered in this investigation, specific Mach number and altitude limits were established and dimensional designs were identified. Utilization of these codes and representative results are documented in this paper. Ó 2005 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail address:
[email protected] (R.J. Hartfield).
0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.003
1540
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1. Introduction The design of a scramjet powered missile is very complex and has been the subject of investigations for over four decades [1–9]. In most cases, designers have tried to take advantage of the vehicleÕs shock structure in order to increase range and/or payload. This leads to a ‘‘wave rider’’ design. It would also be a significant advantage if the design could fly at or near maximum lift to drag ratio throughout the trajectory. If the flight were at constant altitude, the design problem becomes nearly deterministic. That is, with a few assumptions, the vehicle shape that will produce minimum drag and maximum lift may be determined directly. However, the design problems under consideration for many practical applications require extensive variability in altitude and Mach number so that a ‘‘wave rider’’ design is not practical. Instead, an axisymmetric shape is assumed and shape variations are dictated by the internal ballistics of the system. Additionally, it will be required that the scramjet powered missile be launched with an external propulsion source, such as a solid rocket, to some ‘‘take-over’’ Mach number and altitude. At that point, the scramjet engine is started and propels the system until a minimum Mach number and or altitude preclude further operation of the engine. Currently, the success rate for producing thrust in flight conditions using scramjets is unclear. Because of the hypersonic flow regime, it is impossible to separate the external aerodynamics from the internal fluid flow path. Consequently, the system must be designed as a symbiotic collection of interconnected pieces. It is clear, however, that the design must obey the basic laws of conservation of mass, momentum and energy and that the equations of motion applied to the fluid internal to the design must finally result in a positive thrust solution if the system is to experience prolonged powered flight.
2. Basic assumptions Scramjet system models are based on a set of assumptions and, as with all models, the integrity or accuracy of the predictions is no better than the assumptions. For the present case, a combined axisymmetric and onedimensional flow path system is assumed so that cross flows, boundary layer gradients, and the details of the mixing of the air with the fuel are not considered in detail for this preliminary design level optimization. The framework for the modeling has been established in the form of elementary predictive models, which are suitable for demonstrating the optimization process. The working fluid, air, is modeled as a thermally perfect gas so that the ideal gas equation of state holds. Additional assumptions are discussed as follows.
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1541
2.1. Specific heat In addition to the axisymmetric assumption previously mentioned, it is clear that calorically perfect assumptions such as constant specific heat and the concept of constant total temperature are not valid for scramjet systems. A thermally perfect empirical relationship for the determination of specific heat, taken from Mattingly [1], is used here. The equation that defines the specific heat for air is C P ¼ a0 þ a1 T þ a2 T 2 þ a3 T 3 þ a4 T 4 þ a5 T 5 þ a6 T 6 þ a7 T 7 ;
ð1Þ
where the constants are a0 = 0.25020051, a1 = 5.1536879d-5, a2 = 6.5519486d-8, a3 = 6.7178376d-12, a4 = 1.5128259d-14, a5 = 7.6215767d18, a6 = 1.45267720d-21, a7 = 1.0115540d-25. The units for CP are (BTU/Lbm °R). The ratio of specific heats, c, is therefore also a function of temperature and may be expressed as c¼
CP ; CP R
ð2Þ
where R = 0.0685463144 (BTU/Lbm °R). 2.2. Mass conservation In steady fluidic problems, the conservation of mass is implied. However, when the fluid chemistry is changing, care must be taken to ensure that all mass is accounted for. For the present case, this fundamental concept was used as a governing equation and in all cases, all mass was accounted for. 2.3. Conservation of energy Since the specific heat is dependent on temperature, the calorically perfect assumption of a constant total temperature for adiabatic flow problems is not valid. Instead, it is assumed that the total enthalpy is constant. The total temperature is computed using the energy-based definition and varies somewhat along the flow path even for portions of the flow where there is no heat addition. It is assumed that there is not heat transfer for the entire scramjet flow path. The heat transfer rates are not zeros but the adiabatic assumption is reasonable given the flow velocities and residence times. Air enters the inlet at a spatially constant total enthalpy, dependent on the free stream Mach number and altitude. The air then passes through a shock structure and finally enters the combustion chamber. The fuel, hydrogen, is mixed with the air and the total enthalpy of the mixture is increased. This high temperature mixture is assumed to then expand through the nozzle section and exit the system at a lower temperature and pressure but with a total enthalpy of
1542
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
the air–hydrogen mixture. As previously indicated, the specific heat is allowed to change everywhere along the flow path as a function of the local static temperature. 2.4. Momentum The momentum equation is used in the development of the one-dimensional analysis and for the prediction of thrust.
3. Flow modeling The flow field involved in this particular scramjet design is composed of an inlet section, a shock section, a combustor section and an exhaust nozzle. In general, the flow field is an annular shape with the annular ‘‘thickness’’ much less than the local radius. For this reason, the flow field may be treated as two dimensional. 3.1. Basic geometry A schematic of the basic geometry of the hypersonic missile is shown in Figs. 1 and 2, illustrating a design that is axisymmetric. The nose section consists of a single conical nose spike which directs the attached conical shock to the nacelle inlet ‘‘lip’’ when operating at some design altitude and Mach number. One could also use an isentropic ramp for the nose section and some gains would be realized in a reduction in total pressure losses. However, a conical nose structure is a conservative design approach and the utilization of an isentropic ramp would only improve the on design performance. Off design performance for an isentropic ramp is uncertain and the smaller incident cone angle will present additional heat transfer problems for an actual flight vehicle. The method used to generate the solution for a conical shock follows closely that of Ref. [4]. The model assumes a single conical ramp inlet followed by a capture duct or nacelle. The conical spike directs the flow through an annulus to the combustor section. The flow then passes through the combustor and enters an internal nozzle expansion section configured such that the nozzle exit area is equal to the scramjet outer diameter.
Fig. 1. Schematic of a typical scramjet design.
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1543
Fig. 2. Schematic of a scramjet showing the conical shock at the design point.
The inlet is modeled in considerable detail using the conical shock model. The conical shock model and one-dimensional analyses of the bypass duct, combustor and nozzle, are used first to determine the scramjet geometry based on the design point Mach number and altitude. The exit conditions of the inlet can then be calculated for both sub-critical and supercritical flight conditions. 3.2. Cone shock solution The most convenient coordinate system for analyzing an inviscid conical shock is a spherical system centered at the origin of both the physical cone and the shock. The solution methodology followed in this work is outlined in detail in Ref. [4]. The effects of the viscous boundary layer and the angle of attack are neglected in this analysis. The flow is considered to be uniform and supersonic upstream of the conical shock induced by a conical body. The region of the flow between the conical shock and the conical body has a variable velocity field and thermodynamic state. In this region, the velocity in both the r and h directions and the thermodynamic conditions are functions of the angle h only. The flow properties are therefore both axisymmetric and constant along radial rays. The resulting equations as presented in Ref. [4] provide the basis for the solution of the region between the shock and the solid conical surface. The final equation is the well known Taylor MacColl equation which may be written as " 2 # c1 2 dV dV r d2 V r r 2 V max V r cot h þ 2V r þ 2 dh dh dh2 2 dV r dV r dV r d V r Vr þ ¼ 0. ð3Þ dh dh dh dh2 For conical flows, this equation is an ordinary differential equation in h for the variable Vr. Once Vr is determined, Vh can be determined from Eq. (4). Vh ¼
oV r : oh
ð4Þ
1544
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
There is no closed formed analytical solution to this equation and it must be solved numerically. For numerical solution techniques, it is convenient to define a non-dimensional velocity as follows: V0
V : V max
Rewriting Eq. (3) in terms of the non-dimensional velocity yields: " 0 2 # c1 dV r dV 0r d2 V 0r 02 0 1Vr cot h þ 2V r þ 2 dh dh dh2 dV 0r 0 dV 0r dV 0r d2 V 0r Vr þ ¼ 0. dh dh dh dh2
ð5Þ
ð6Þ
The numerical procedure adopted for solving this equation for a given cone angle and flight Mach number involves solving the inverse problem iteratively. To solve the inverse problem, the following methodology is adopted: 1. A cone shock angle, hs, is assumed for a given free stream Mach number, M1. From these two parameters, the flow (deflection) angle and Mach number immediately downstream of the oblique shock can be obtained using the oblique relations for planar shocks. The relationship for the flow angle is M 2 sin2 hs 1 tan h ¼ 2 cot hs 2 1 ; ð7Þ M 1 ðc þ cos 2hs Þ þ 2 where h is the flow angle immediately downstream of the shock and c is the ratio of specific heats. Once the flow angle immediately down stream of the shock is determined, the downstream Mach number can be determined from M2 ¼
M n2 ; sin ðhs hÞ
ð8Þ
where M n2 is the downstream Mach number determined from the normal shock relations using an upstream Mach number for a flow direction perpendicular to the shock. 2. From the Mach number and flow angle immediately down stream of the shock, non-dimensional radial and angular components of flow velocity, V 0r and V 0h , can be determined for the position immediately down stream of the shock. Note that V 0 is strictly a function of Mach number and can be determined from 12 V 2 0 ¼V ¼ þ1 : ð9Þ V max ðc 1ÞM 2
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1545
3. Using the values of V 0r and V 0h calculated in step 2 as boundary conditions, Eq. (6) can be solved numerically in h, marching away from the shock. For this work, a fourth-order Runge–Kutta integration routine is used to solve Eq. (6) for the values of V 0r at each Dh, recognizing that V 0h is the derivative of V 0r as follows from Eq. (4). It is very important to recognize that, for the problem as formulated above, V 0h will take on a substantially negative value at the boundary and will steadily increase with h (its magnitude will decrease). 4. As the integration progresses, the value of V 0h will pass through 0.0 for some value of h. This value of h corresponds to hc, the physical cone angle. At this location, the velocity V 0 is equal to V 0r and the Mach number near the surface of the cone can be determined from Eq. (9). The velocity field has been determined in the integration process and the pressure and temperature can be determined at each point from the isentropic flow relations since the only entropy change occurs across the oblique shock. 5. The cone angle determined in step 4 is not likely to be the cone angle selected for the missile. Hence, steps 1–4 are followed in an iterative fashion until the cone angle calculated in step 4 corresponds to the cone angle of the inlet. The above discussion addresses the design point operation. For a single case simulation, the design point geometry is determined as outlined above. As the engine is flown, the inlet will operate in either the sub-critical, critical or super critical mode. The analysis for these modes is discussed in a later section. 3.3. Conical shock solution verification To gain some confidence in the cone shock solution, the results were validated by using tabulated values for the cone shock solution found in NACA Publication 1135 for a specific heat ratio of 1.4. The comparison between the calculated values and the NACA results were identical for the cases presented in 1135. 3.4. Inlet analysis Once the conical flow solution is determined, the properties of the air behind the conical shock are found from standard oblique shock relations. The initial conical shock is either a ‘‘design point’’ case or an off-design case. In either case, the flow must be decelerated still further before entering the combustor. This deceleration is accomplished by shocking the flow down through a series of internal oblique shocks. The number of internal shocks can range from one to several but one can determine quickly that a single oblique shock creates too much total pressure loss and consequently, two or more internal shocks are required. For the present case, two shocks were considered.
1546
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
3.4.1. Design point operation When the inlet is operating at the design point, the conical shock impinges on the nacelle leading edge as shown in Fig. 3, and the capture area for the inlet is the frontal area of the nacelle. Behind the conical shock, the air remains parallel to the cone at the conical surface and parallel to the shock immediately downstream of the shock. The inlet mass flow rate is simply dependent on the capture area outlined by the interior dimension of the inlet. That is, for the design case, rffiffiffiffiffiffiffiffi cg m_ a ¼ pa Acap M a ; ð10Þ RT a where pa and Ta are the free stream static pressure and temperature, Ma is the flight Mach number, R is the gas constant for air (R = 53.34 ft/°R), Acap is the capture area defined by Acap ¼ pr2inlet . Immediately downstream of the shock, the Mach number is not a constant and thus computation of the mass flow rate is dependent on the appropriate ‘‘local’’ conditions. In any case, the mass flow rate in this region (region 2 in Fig. 6) may be written as rffiffiffiffiffiffiffiffi cg m_ 2 ¼ p2 Arec M 2 : ð11Þ RT 2 For the conical shock, the Mach number immediately behind the shock is used to determine the fluid properties in this region. The flow area, Arec, is determined from the ‘‘radial’’ geometry as Acap Arec ¼ 2 ðcos hN cos bÞ; ð12Þ sin2 b where b is the conical shock angle and hN is the half angle of the cone.
Nacelle Vθ ck Sh o
M = MDesign
Vr ε
o se β N θc
C one
Spike
Fig. 3. Design point for a conical shock inlet.
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1547
Eq. (11) is then used to determine the mass flow rate in region 2. Since the Mach number used in Eq. (11) is known to be incorrect, it is adjusted until the mass flow rate from Eqs. (11) and (10) agree. This is equivalent to finding an average Mach number in this region and using this average value to determine the fluid properties in region 2. 3.4.2. Sub-critical operation A schematic of sub-critical operation of a conical inlet is shown in Fig. 4. For sub-critical operation, the conical shock extends out in front of the nacelle lip. The mass flow rate of the air through the engine is decreased and may be determined using the ratio of the annular areas. If the Mach number at the inlet is assumed to be constant, then the ideal mass flow rate is determined as in the previous design case, so that rffiffiffiffiffiffiffiffi cg m_ aideal ¼ pa Acap M a : ð13Þ RT a The total area associated with this ideal flow rate is a segment of the surface of a ‘‘sphere’’ with its origin at the apex of the cone. That is, the flow field appears as if it emanates from the conical apex with a radial speed Vr. This total area is found from Eq. (12), where b is the actual shock angle for the sub-critical case. The actual flow area is ! Acap ð14Þ Aflow ¼ 2 cos hN cos hLip : 2 sin hLip Thus the actual flow rate is m_ offdesign ¼ m_ aideal
Atotal : Aflow
ð15Þ
Nacelle
k
oc Sh
M = MSub-critical
Vθ
Vr
ε
ose β N θc
Cone
Spike
Fig. 4. Sub-critical operation of a conical inlet.
1548
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
Once, the mass flow rate is determined and the inlet parameters have been established behind the shock, the solution proceeds in exactly the same manner as for the design case. 3.4.3. Supercritical operation A schematic of supercritical operation of a conical inlet is shown in Fig. 5 and for this supercritical operation, the conical shock extends inside the front edge of the nacelle lip. The capture area for this mode of operation is the same as design point capture area and the nozzle mass flow rate can be calculated as it was in the design point case. There are additional shocks that should be considered but are omitted for the present analysis. In lieu of these added shocks, the shock angles inside the ‘‘shock section’’ are altered to satisfy continuity. This approach will result in small changes in the actual total pressure losses, but is a good approximation to a ‘‘fixed’’ geometry case. As shown in Fig. 5, there is an additional external shock attached to the nacelle lip which will cause an additional drag increment. In the present model, this increment in drag is neglected. This completes the conical inlet section for both the design case and the sub-critical and supercritical cases. After the flow passes through the conical section, its speed must be further reduced in the ‘‘shock’’ section. 3.5. Internal shock section The internal shock section is designed to slow the flow down before it enters the combustion section. For the present design cases, a two shock system was used. A schematic of this section is shown in Fig. 6. This completes the conical inlet section for both the design case and the sub-critical and supercritical cases. Nacelle
Vθ k hoc
S
β
M = MSuper-critical
Nose
Vr Cone
Spike
θc
Fig. 5. Supercritical operation of a conical inlet.
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1549
Nacelle Vθ Vr
Shoc k -1 3 2
ε
Nose
Cone
Shock - 2
4
Spike
Fig. 6. Schematic of internal shock section.
The shock separating Sections 2 and 3 and the shock separating Sections 3 and 4 (see Fig. 6) are assumed to be planar. Actually, these shocks are curved sections of an annular segment, but because the annular radius is much smaller than the average radius, a planar shock assumption is a good approximation. The deflection angles, which would be fixed in the actual hardware, are assumed to be ‘‘variable’’ in the present case so that continuity and energy can be satisfied. That is, even after the design point has been established, the deflection angles creating these shocks are assumed to vary in order to satisfy continuity and energy. For the design case, the deflection angles for these two shocks are chosen by the GA routine. It is not intuitive to the designer what the value of the first angle should be relative to the total turn angle. It is initially assumed that the total turn angle for both these shocks is the same as the nose cone angle. That is, if the flow at the entrance to this section were parallel to the conical spike, then to turn the flow back parallel to the free stream would require a total turn angle of hN. The flow is, however, not quite parallel to the nose cone and thus the total turn angle is not quite hN. Hence, the assumption of allowing the angles to vary in order to satisfy continuity and energy is a valid approximation. Later computations did, in fact, confirm that the total turn angle is very close to hN. The determination of the turn angles to satisfy continuity and energy is an iterative procedure for which there is no direct computation. First, a value for the first deflection angle is assumed and the ‘‘planar’’ shock properties are determined for the fluid properties in Section 3. The second shock angle is then assumed and the properties are determined in Section 4. At this point, the mass flow rate is found since the annular geometry may be determined from the wall
1550
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
deflection angles. If this is not the same as the mass flow rate at Section 2, then new wall deflection angles must be assumed. An iterative technique was developed to solve this problem and solutions were found for cases in which a solution existed. Once shock angles have been found to satisfy continuity and energy, the computation moves to the combustor section. Several checks are made at the entrance to the combustor to insure that the fluid possesses sufficient Mach number properties and low enough static temperatures to allow combustion to take place. For the scramjet case, this means that the Mach number must be around 3.0 or less with a static temperature around 1000 °R or less. The exact values of these variables have been determined by the physics of the problem and the design variables initially chosen by the GA at the start of the design. 3.6. Combustor model The scramjet combustor is modeled in a one-dimensional sense. A fuel heating value, fuel density, combustion efficiency, and the fluidic conditions at the exit of the shock section are used as inputs to the combustor. For this effort, hydrogen is used as the fuel and the combustion efficiency is assumed. The fuel to air ratio for a hydrogen–air mixture is determined from the ideal combustion of hydrogen and air. The chemical equation for the stoichiometric combustion process is 2H2 þ ðO2 þ 3:76N2 Þ ! 3:76N2 þ 2H2 O
ð16Þ
and the mass ratio for stoichiometric conditions is far ¼
MassH2 2MwH2 ¼ ¼ 0:0145687: Massair MwO2 þ 3:76MwN2
ð17Þ
It is further assumed that, for the design case, the ideal fuel–air mixture ratio is used and for all other off-design cases, the fuel–air ratio is controlled by a fuel flow rate controller to dispense fuel as necessary to maintain the desired operating conditions. 3.6.1. Combustor length model The required length of the combustor to fully combine the hydrogen fuel with the oncoming air stream is an area of uncertainty. No attempt is made, in the present case, to define a detailed combustor length model. Instead, an empirical equation is used to produce a reasonable combustor length that is based in part on the combustors used in ramjet designs. The model presented by Mattingly [1] was modified and utilized for the higher speed jets in a scramjet engine. It is recognized that a more accurate combustor length would be a desirable feature and should be pursued in future modeling endeavors. For the present case, the equation defining the combustor length is
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
LCombustor ¼ 0:5M 24
0:714 68843:0P 04 pffiffiffiffiffiffiffi : T 05
1551
ð18Þ
Even though T05 is not known with precision, it may be estimated to within a few percent. In any case, the combustor length is weakly dependent on the exit temperature and more strongly dependent on the entrance Mach number. The constants in the equation were chosen from satisfactory performance of ramjet engines already in existence. This equation assumes that the temperature is in °R and the pressure is in psia. A plot of Eq. (18) is shown in Fig. 7, illustrating the dependency of the combustor length on incoming pressure and Mach number. For this particular plot, it was assumed that the static temperature at the exit of the combustor was 3800° R and that the average value of c for this section was 1.3. 3.6.2. Combustor fluidic model Inside the combustor, the rise in enthalpy is assumed to follow a curve in which the temperature initially increases rapidly followed by a gradual increase to its final value. The empirical equation used to describe the total temperature in the combustor is 3 2 x tan1 20:0 LComb 5: T 0 ¼ T 04 þ ðT 05 T 04 Þ4 ð19Þ tan1 ð20:0Þ In this equation, the x variable is 0.0 at the entrance and is equal to LComb at the end of the combustor. This equation is empirical but describes a rapid total temperature buildup which emulates the combustion of H2 in air.
300 Po4=50
L_combustor
250 200 Po4=100 150
Po4=150 Po4=200
100
Po4=300 Po4=500
50
Po4=1000
0 1
2
3 Mach No. (4)
Fig. 7. Combustor length model.
4
5
1552
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
The differential equations which describe the flow field inside the combustor include the effects of (1) heat addition (Raleigh flow), (2) friction (Fanno flow), (3) area change, and (4) changes in specific heats. The fundamental equation describing these changes for a one-dimensional analysis is written in terms of the Mach number and is expressed as M 2 1 dA M 2 1 þ cM 2 1 þ c1 M2 1 dM 2 2M 2 1 þ c1 2 2 þ ¼ A dx T0 dx ð1 M 2 Þ ð1 M 2 Þ c1 4 2 2 dT 0 4f cM 1 þ 2 M 1 M dc : ð20Þ þ D dx c dx ð1 M 2 Þ In this equation, f is the skin friction coefficient and D is the hydraulic diameter. It is the solution of this differential equation that governs the properties of the flow inside the combustor. The equation describing the cross-sectional area of the combustor is assumed to be quadratic in nature, so that the area change term is 1 dA A4 k 1 x ¼ 1 ; ð21Þ A dx AðLC Þ LC where k1 is a ‘‘variable’’ chosen by the GA routine. The total temperature change term is found from the derivative of Eq. (19) which may be written as 1 dT 0 ðT 05 T 04 Þ 20LC 1 ¼ ð22Þ 1 2 2 T 0 dx tan ð20:0Þ LC þ 400x T 0 and the change in c is estimated to be 1 dc c3800 cT 4 ¼ : c dx cLC
ð23Þ
The total temperature is related to the total enthalpy so that h0 = CPT0, where CP is a function of local static temperature defined by Eq. (1). Eq. (20) is a first order equation in M2 and is solved using a fourth order Runge–Kutta integration scheme. The result is a definition of the Mach number at each station along the combustor. From the Mach number at each point, pressure, static temperature, enthalpy and density are available through the following equations. The fuel mass flow rate is m_ f ¼ m_ a ðfar Þ ðlbm =sÞ; where far is the fuel to air ratio so that the total mass flow rate becomes m_ T ¼ m_ a þ m_ f ðlbm =sÞ:
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1553
5
Fuel
e
Fig. 8. Schematic of the nozzle section.
The static pressure is sffiffiffiffiffiffiffi m_ a ð1 þ far Þ RT p¼ : AM cgc
ð24Þ
The static temperature, T, and the total pressure, P0, are found from standard isentropic relations. For the density, the ideal gas equation of state is used. At the end of the combustor, the final flow conditions are determined and a check is made on the static temperature. The maximum allowable static temperature at any point in the combustor is 3800 °R. Above this value, the gas will begin to dissociate and energy is pumped into non-translational modes of energy which cannot be recovered to produce thrust. In addition to the fluid properties, the cross-sectional flow area is found and the final expansion through the nozzle begins. 3.7. Nozzle section In the nozzle, the flow expands to the exit velocity at the physical exit plane. For the present design case, it is assumed that the flow expands around an internal plug as shown in Fig. 8. That is, there is no external nozzle extending outside the nacelle. The initial curvature is a segment of a circle and the concave section is a segment of a parabola. It is assumed that no shocks exist in the nozzle section. 4. Computation of thrust The thrust is determined from momentum considerations from the inlet conditions to the conditions at the exit plane. The equation for this term in the thrust equation is T hr1 ¼ ðm_ e me m_ a ma Þ=gc : ð25Þ
1554
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
This is usually the largest term in the thrust equation, but other terms also make significant contributions. The pressure term is T hr2 ¼ ðpe pa ÞAe :
ð26Þ
These are the two basic terms in the thrust equation and other terms contribute only to the overall drag. It is assumed that the drag terms consist of a pressure drag term and a skin friction term. The pressure term results from the variation of the surface pressure on the nacelle surface. Since the exterior shape of the nacelle section is a section of a circular arc, the pressure may be determined exactly if an isentropic expansion is assumed over the exterior section of the nacelle. This pressure drag is determined by integration of the pressure distribution on the nacelle surface in the axial direction. The skin friction is determined by computing the local skin friction coefficient along the surface (function of the local Reynolds number) and summing the resulting forces from the nacelle to the end of the missile. For these computations, it was assumed that the boundary layer was turbulent. Even though the boundary layer may, in fact, be laminar at high altitudes, a turbulent model is a conservative approach. Finally the net thrust is Thrust ¼ T hr1 þ T hr2 Dpress Dfric ;
ð27Þ
which includes the dominant thrust and drag terms in the system. The only term missing is the frictional losses on the entrance cone.
5. Results In this analysis, the scramjet system is designed a Mach number for which the conical nose shock impinges exactly on the lip of the nacelle. This flight condition and resulting geometry is termed the design point. Once the design point has been established, the off-design performance is determined at various altitudes and Mach numbers in order to clearly identify the flight envelope for that particular case. For the present analysis, a design altitude and design Mach number were chosen and the design geometry for that case was determined. This geometry was determined with using a genetic algorithm with maximum specific impulse as the goal. The specific impulse is defined as I SP ¼
Thrust lbsF : ... lbsm =s m_ Fuel
ð28Þ
The variables used by the GA to maximize the specific impulse are (a) the conical nose angle (b) the two internal shock angles, (c) the area change variation in the combustor and (d) the needed fuel to air ratio. Consequently, for a chosen design altitude and Mach number and the determination of these four
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1555
Table 1 Run matrix for scramjet designs Case
Md
Altd
hN
DInlet
Q, psia
I SPMax
ThMax
1 2 3 4 5 6 7 8
7 9 7 7 7 5 5 10
80K 100K 100K 80K 100K 66K 100K 110K
10.97 8.84 10.26 11.68 10.97 15.22 15.22 6.71
12 12 12 20 8 12 12 12
1974 1279 774 1974 774 1972 395 1011
2086 1534 1434 2360 995 2767 2994 1368
438 161 117 1407 36 831 167 100
variables, the specific impulse is maximized. A summary of the cases considered are presented in Table 1. The resulting specific impulse and thrust for these design cases are also presented. For each of these designs, the off-design performance was determined by considering altitudes ranging from 40,000 ft to 118,000 ft and Mach numbers at each altitude ranging from M = 5.0 to M = 11.0. At each Mach number—altitude combination, the specific impulse and net thrust were determined. Restrictions or constraints were placed on the resulting solutions by limiting the maximum static temperature anywhere in the engine to 3800 °R. The maximum dynamic pressure was set at 2000 lb/ft2 and the Mach number was not allowed to go below M = 1.1 anywhere in the system. If any of these conditions were not met, the run was terminated and another Mach number— altitude combination was chosen. Implicit in these computations was the requirement that the overall thrust be positive and positive thrust resulting from the momentum term in the thrust equation. One note of clarification is in order. If a design point is optimized to produce maximum specific impulse, this does not mean that higher specific impulses are not possible for off-design cases. That is, a high Mach number design point can indeed lead to higher specific impulse performance at lower Mach numbers in an off-design flight condition but would remain lower than the specific impulse of a scramjet optimized for the lower Mach number. This condition is underscored in this analysis by the fact that the geometry was held constant at the design point. The total enthalpy was forced to remain constant from the inlet to the beginning of the combustor. This means that the two shock angles were allowed to vary somewhat even for off-design cases in order to satisfy this requirement. This is equivalent to allowing the fluid to pass through multiple internal shocks (instead of just two) in order to satisfy continuity and energy as would actually occur in a fixed geometry engine. The resulting flight envelopes are presented in Figs. 9–24 for the cases as presented in Table 1. It is clear from these plots that for any given design point, the operating envelope is very narrow. That is, the combination of Mach
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 3000 Dinlet = 12 in Altd = 80K ft Md = 7.0 N = 10.97 deg far = 0.01449
2500 82K ft
Design Point
2000
Isp
86K ft
1500
94K ft 102K ft
1000 500 0 6
6.5
7
7.5
8
Mach No.
Fig. 9. Case 1, Mach no. vs ISP, Md = 7.0, Altd = 80K ft.
500 Design Point
Dinlet = 12 in Altd = 80K ft Md = 7.0 N = 10.97 deg far = 0.01449
82K ft
Thrust
400
86K ft
300
200 94K ft 100 102K ft 0 6
6.5
7
8
7.5
Mach No.
Fig. 10. Case 1, Mach no. vs thrust, Md = 7.0, Altd = 80K ft.
3500 84k ft
Dinlet = 12 in Altd = 100K ft Md = 9.0 N = 8.84 deg far = 0.01449
90k ft
3000
96k ft 2500 102 k ft 2000
Design Point
Isp
1556
1500 118 k ft 1000 500 0 7
7.5
8
8.5
9
9.5
Mach No.
Fig. 11. Case 2, Mach no. vs ISP, Md = 9.0, Altd = 100K ft.
10
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 700 Dinlet = 12 in Altd = 100K ft Md = 9.0 N = 8.84 deg
600 84K ft
500
far = 0.01449
Thrust
90K ft
400 96K ft
300
102K ft
Design Point
200 100
118K ft
0 7
7.5
8
8.5
9
9.5
10
Mach No.
Fig. 12. Case 2, Mach no. vs thrust, Md = 9.0, Altd = 100K ft.
3000
Dinlet = 12 in Altd = 100K ft Md = 7.0 N = 10.26 deg far =0.0146
72K ft 90K ft
2500
118K ft
Isp
2000 1500
Design Point 1000 100K ft
500 0 4
5
6
7
8
Mach No.
Fig. 13. Case 3, Mach no. vs ISP, Md = 7.0, Altd = 100K ft.
800 Dinlet = 12 in Altd = 100K ft Md = 7.0 N = 10.26 deg far = 0.0146
700 72K ft
600 Thrust
500 400 84K ft
300 90K ft
200 100K ft
100
Design Point
118K ft
0 4
5
6
7
Mach No.
Fig. 14. Case 3, Mach no. vs thrust, Md = 7.0, Altd = 100K ft.
8
1557
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 2500 Dinlet = 20.0 in Altd = 80K ft Md = 7.0 N = 11.68 deg far = 0.0146
Design Point 2000 84 k ft
ISP
1500
1000
500 118 k ft
92 k ft 106k ft
0 6
6.5
7
7.5
8
8.5
9
Mach No.
Fig. 15. Case 4, Mach no. vs ISP, Md = 7.0, Altd = 80K ft.
1600 Design Point
Dinlet = 20 in Altd = 80 K ft Md = 7.0 N = 11.68 deg far = 0.0146
1400
Thrust
1200
84k ft
1000 800 600 92k ft
400 106k ft
200 118k ft
0 6
6.5
7
7.5
8
8.5
9
Mach No.
Fig. 16. Case 4, Mach no. vs thrust, Md = 7.0, Altd = 80K ft.
3000
Dinlet = 8.0 in Altd = 100 K ft Md = 7.0 N = 10.97 deg far = 0.0146
2500 2000
Isp
1558
1500 118k ft
82 k ft 96 k ft
1000 Design Point 500 0 4
5
6
7
Mach No.
Fig. 17. Case 5, Mach no. vs ISP, Md = 7.0, Altd = 100K ft.
8
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 300
Dinlet = 8.0 in Altd = 100 K ft Md = 7.0 N = 10.97 deg far = 0.0146
72k ft
250
Thrust
200 150
82k ft
100 96k ft
50 118k ft
Design Point
0 4
5
6
7
8
Mach No.
Fig. 18. Case 5, Mach no. vs thrust, Md = 7.0, Ald = 100K ft.
3000
Design Point 100K ft
2000
Isp
Dinlet = 12 in Altd = 66 K ft Md = 5.0 N = 15.22 deg far = 0.0146
70K ft
86K ft
2500
118K ft
1500 1000 500 0 4
5
6
7
8
Mach No.
Fig. 19. Case 6, Mach no. vs ISP, Md = 5.0, Altd = 66K ft. 900
Design Point
Dinlet = 12 in Altd = 66 K ft Md = 5.0 N = 15.22 deg far = .0146
800 70K ft
700
Thrust
600 78K ft
500 400
86K ft
300 200
100K ft
100
118K ft
0 4
5
6
7
Mach No.
Fig. 20. Case 6, Mach no. vs thrust, Md = 5.0, Altd = 66K ft.
8
1559
1560
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 1800 Dinlet = 12 in Altd = 100 K ft Md = 5.0 eN = 15.22 deg far = .0146
1600 70 k ft
Isp
1400
Design Point
1200 116 k ft
1000 800
104 k ft
600 4
4.5
5 Mach No.
5.5
6
Fig. 21. Case 7, Mach no. vs ISP, Md = 5.0, Altd = 100K ft.
200 Dinlet = 12 in Altd = 100 K ft Md = 5.0 eN = 15.22 deg far = .0146
180 160
Thrust
140 120 100
70 k ft
80 Design Point
60
40
86 k ft
20
116 k ft
104 k ft
0 4
4.5
5
5.5
6
Mach No.
Fig. 22. Case 7, Mach no. vs thrust, Md = 5.0, Ald = 100K ft.
numbers and altitudes at which positive thrust is obtained is very limited. On the lower Mach number end, the pressures become very large and structural limits are exceeded. At the upper Mach number end, the internal temperatures become very large. At low altitudes, the pressure and density are high and at high altitudes, the mass flow rates through the engine become very small and the thrust that is generated is very small. The performance of the engine is dependent on the square of the Mach number and the atmospheric pressure is exponentially dependent on altitude. Consequently, it is not surprising that the operating envelope is very narrow. The size of the operating envelope is highly dependent on the design point. In Figs. 25 and 26, the operating envelope is plotted as a function of Mach number and altitude. Some envelopes are quite small as in Case 6 and this case
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 3500 90 k ft
Dinlet = 12 in Altd = 110 K ft Md = 10.0 N = 15.22 deg far = 0.0146
3000 100 k ft 108 k ft
2500
Isp
118 k ft
2000
1500 Design Point
1000
500 8
8.5
9
9.5
10
10.5
11
11.5
12
Mach No.
Fig. 23. Case 8, Mach no. vs ISP, Md = 10.0, Altd = 100K ft.
600 Dinlet = 12 in Altd = 100 K ft Md = 10.0 N = 15.22 deg far = .0146
90 k ft
500
Thrust
400 94 k ft
300 200 Design Point
100 100 k ft
118k ft
0 8
8.5
9
9.5
10
10.5
11
11.5
12
Mach No.
Fig. 24. Case 8, Mach no. vs thrust, Md = 10.0, Altd = 100K ft.
130 120 110
Alt
100
Case 6
90
Case 1
Case 8
80 70 60 50 40 4
5
6
7
8
9
10
11
Mach No.
Fig. 25. Flight envelope for selected design point cases.
12
1561
1562
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563 130 120 110
Alt
100
Case 7
90
Case 2 Case 3
80 70 60 50 40 4
5
6
7
8
9
10
11
12
Mach No.
Fig. 26. Flight envelope for selected design point cases, Altd = 100K ft.
would result in a very limited operating Mach number range. A more attractive case is either Case 2 or Case 3 where the flight envelope is relatively large. Even in this case the take-over Mach number is very large and the maximum attainable Mach number is less than 10 even for Case 2. The conclusion that can be drawn from these figures is that additional design variables will be required if a relative low take-over Mach number is required along with a high altitude, high Mach number cruise condition.
6. Conclusions and recommendations The work presented in this paper represents a fundamental preliminary design strategy based on sound, first order modeling principles and it represents the first step taken toward a global optimization approach to scramjet design. In this design effort, the continuity and energy equations are satisfied and the momentum equation is used to predict thrust. A variable specific heat, essential to an analysis of this nature is included as well as geometric variations required to define engine operation. This demonstration has shown that the ‘‘neighborhood’’ of the design space which is likely to produce workable designs for a given requirement can be identified using a genetic algorithm connected to a suitable, yet practical preliminary design tool. For a device such as a scramjet, which will inherently operate on the edge of physical viability, design optimization beginning with an approach such as this may very well represent the only viable approach to producing an operating system. Several straightforward improvements and expansions could be made to the preliminary design tool without jeopardizing its usefulness as an efficient performance predictor. Such an effort is likely to represent the next step in the scramjet research program begun with this effort. These improvements include the addition of a boundary layer model, improvement to the combustor length model, the inclusion of multiple shocks in the ‘‘shock’’ section, the inclusion of
R.J. Hartfield et al. / Appl. Math. Comput. 174 (2006) 1539–1563
1563
some external heating modeling, and the inclusion of some external aerodynamics modeling. Additionally, the calibration and some more detailed validation of the preliminary design tool using CFD would improve the quality of the results in many cases. Additionally, the inclusion of variable geometries for the inlet, combustor and nozzle would improve the versatility of the optimization effort.
Acknowledgments The authors wish to thank the US Army Aviation and Missile Command at Redstone Arsenal, AL for their help and support of this effort.
References [1] [2] [3] [4] [5] [6] [7] [8] [9]
J.D. Mattingly, Elements of Gas Turbine Propulsion, McGraw Hill, 1996, pp. 820–822. D. Kuchemann, The Aerodynamic Design of Aircraft, Pergamon Press, Oxford, 1978. J.D. Anderson, Hypersonic and High Temperature Gas Dynamics, McGraw-Hill, 1989. J.D. Anderson, Modern compressible flow with historical perspective, Series in Aeronautical and Aerospace Engineering, McGraw-Hill, 1990. C. Park, Nonequilibrium Hypersonic Aerothermodynamics, John Wiley, 1990. W.H. Heiser, Hypersonic Airbreathing Propulsion, AIAA Press, 1994. J.J. Bertin, L.J. Smith, Aerodynamics for Engineers, third ed., Prentice Hall International, 1998. E.T. Curran, S.N.B. Murthy (Eds.), Scramjet Propulsion, AIAA Progress in Astronautics and Aeronautics, vol. 189, 2000. L.H. Schindel, Factors affecting performance of Ramjet and Scramjet powered vehicles, NSWC TR 87-190, Naval Surface Weapons Center, Strategic Systems Department, Dahlgren, VA, May 1987.