Acta Astronautica 61 (2007) 866 – 872 www.elsevier.com/locate/actaastro
The optimization of rough surface supersonic nozzles D. Vlassova , J.V.C. Vargasa,∗ , J.C. Ordonezb a Departamento de Engenharia Mecânica, Universidade Federal do Paraná, CP 19011, Curitiba, PR, 81531-990, Brazil b Department of Mechanical Engineering and Center for Advanced Power Systems, Florida State University, Tallahassee, FL 32310, USA
Received 11 April 2005; accepted 4 January 2007 Available online 12 April 2007
Abstract In this work, a theoretical optimization study of compressible gases adiabatic flows was conducted for conical nozzles with internal rough surface. A mathematical model was developed for allowing the determination of thrust accounting for nozzle friction and divergence losses. In the first part of the study, it is shown that in rough surface nozzles there are critical values of the nozzle opening half angle, under which, additional increase in gas velocity is not obtained. It is also observed that the minimum nozzle opening half angle value increases as the flow velocity and the nozzle internal surface roughness increase. In the second part of the study, it is shown that there are optimal values of nozzle opening half angle for minimum hydraulic loss across the nozzle, for each combination of the hydraulic resistance caused by the internal surface roughness and the nozzle outlet cross sectional area. © 2007 Elsevier Ltd. All rights reserved. Keywords: Gas dynamics; Supersonic nozzles; Hydraulic resistance; Nozzle outlet cross sectional area
1. Introduction The continuous development of thermal engines that utilize supersonic nozzles to create thrust demand the reduction of nozzle losses to a minimum level. The high gas flow velocity at the divergent part of the nozzle and the roughness of the nozzle internal surface increase considerably friction losses. It is well known that a rocket engine develops maximum thrust when the pressure at the Laval’s nozzle outlet is equal to the atmospheric pressure. As the altitude varies with the rocket trajectory, the atmospheric pressure varies and the engine loses part of its thrust. To compensate such thrust loss it is necessary to increase (or reduce) the length and the area at the ∗ Corresponding author. Tel.: +55 41 3361 3307; fax: +55 41 3361 3129. E-mail address:
[email protected] (J.V.C. Vargas).
0094-5765/$ - see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2007.01.068
nozzle outlet, according to the trajectory, i.e., according to the actual ambient conditions around the rocket position. Therefore, it is necessary to use nozzles that allow a dynamic variation of their geometry. One of the designs that allows the length and the area at the nozzle outlet to vary consists of two parts, a short and a longer one which is the continuation of the former. Such nozzles are manufactured with a flexible material, which is obviously rough. The roughness of the nozzle internal surface, determined by the structure of the material, increases the hydraulic loss even more with the gas expansion increase in the nozzle [1]. In conventional nozzles, in actual operating conditions, the increase in internal surface roughness can be caused by several reasons. In solid fuel rocket engines, the condensed fluid in the liquid phase sticks on the nozzle surface, and changes to the solid phase. The condensed phase particles accumulate on the surface, turning it similar to a sandpaper surface.
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
Another situation that could cause roughness augmentation of the supersonic nozzle is the utilization of a thermal protection ablation polymer to manufacture the nozzle. The nozzle internal surface under the high temperature gas flows is subjected to erosion that augments roughness. In other cases, when special materials are used, the surface roughness is conditioned by the structure of the material itself. It must be noted that the listing of the presented roughness augmentation factors is not complete. In practice, it is likely that a combination of the previously cited factors is responsible for considerably increasing surface roughness, and therefore nozzle friction losses. Conical nozzles have been the subject of research for many years, and many analysis exist in the literature. A few examples related to the present work are herein discussed. Hoffman [2] studied conical nozzles with nonisentropic flow, without considering a variable geometry. Hurst et al. [3] developed experiments and calculated numerically heat transfer coefficients along the hot gas nozzle smooth wall for high-speed boundary layer nozzle flows under engine Reynolds number conditions, presenting a comparison between the experimental and numerical results, showing that the simplification of a local average density in the numerical model leads to good results for transonic and low supersonic flows. Zaccaria et al. [4] performed experiments to determine the three-dimensional flowfield in the nozzle of a low-speed, single-stage, axial-flow turbine, comparing the measurements to the numerical results of an existing Navier–Stokes code with good agreement, basically seeking a better understanding of the nozzle flowfield. Local surface heating was implemented experimentally in a nozzle, by Demetriades [5], showing that increasing surface temperatures decreases the instability downstream, suggesting that the surface heating is important in attenuating nozzle boundary layer instabilities and delaying transition, but not considering surface roughness variation. Another study [6] investigated the possibility of designing an innovative propulsion system for a space vehicle by means of a laser beam based on the repetitively pulsed method, assessing feasibility using an engine with a parabolic nozzle. A mathematical simulation was applied to analyze the effect of the nozzle wall roughness to determine the sensitivity of laser intensity distribution within the focal region to nozzle surface quality, showing that a suitable engine can be designed. A design theoretical analysis and experimental study of supersonic flow of ideal gas around the jet nozzle and flight vehicle with relief surfaces was carried out by Semyonov [7], comparing the measured wave drag data with the calculated
867
theoretical results, with good agreement, but no variation in surface roughness was investigated. Another direction has been pursued in nozzle design to improve performance, which is the use of plug-nozzles, and appropriately designing the plug contour. Following this path, Ito et al. [8] used the method of characteristics to design the plug contour, showing that the altitude influenced base pressure is important to maintain the thrust performance, assuming the chamber pressure is constant during ascent, and that the contoured plug nozzle thrust performance is 5–6% higher than the conical plug nozzle. Therefore, the review of the literature shows that not too much attention has been given to the increase in internal surface roughness in nozzles during operation and its effect on thrust performance. Nozzle optimization has been the objective of scientific research for long. Byington and Hoffman [9] analyzed the base pressure effects on conical thrust nozzle optimization, neither considering a variable geometry nor roughness variation. Addy and White [10] performed the optimization of drag minimums without allowing for a dynamic variation of nozzle geometry. More recent nozzle optimization studies considered supersonic species transport, flow parameters, different nozzle geometries [11], multiple objective design optimization of internal aerodynamic shape operating in transonic flow with genetic algorithms [12], and length optimization of expansion deflection nozzles investigating the effect of various throat parameters on nozzle contour [13]. However, within the knowledge of the authors, the dynamic variation of nozzle geometry and the increase in internal surface roughness have not been investigated simultaneously in a nozzle optimization study previously. In order to provide new information on the aspects discussed in the previous paragraphs, for practical engineering use, the main objectives of the present paper are: (i) to study the variation of the nozzle opening half angle, and its effect on gas velocity, considering the variation of the nozzle internal surface roughness, and (ii) to optimize the values of nozzle opening half angle for minimum hydraulic loss across the nozzle, considering the hydraulic resistance caused by the variation of both internal surface roughness and nozzle outlet cross sectional area. 2. Theory1 One of the most important characteristics of a rocket engine is the thrust developed in vacuum [14]. The 1 Notation is given in Appendix at end of paper.
868
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
absolute value of such thrust does not depend on the atmospheric pressure, and it is determined only by internal parameters of the combustion chamber and the nozzle. Fig. 1 shows a general diagram of the axisymmetric nozzle treated in this study. Assuming unidirectional flow (jet direction), the thrust created in vacuum by the rocket engine is determined by the following equation [14]: P = pc F ∗
2 k+1
1/(k−1) z(a ),
(1)
where P is the vacuum thrust, N; pc is the combustion chamber pressure, Pa; k is the adiabatic ratio, i.e., the ratio of specific heats at constant pressure and constant volume of the gas, cp /cv ; z()=+1/ is a gas dynamics function; = V /V ∗ is the dimensionless velocity; √ V is the gas velocity, m/s; V ∗ = 2(k/(k + 1))RT c is the critical velocity, m/s; a is the nozzle outlet velocity coefficient; R is the gas constant, J/(kg K); Tc is the gas stagnation temperature in the combustion chamber, K; F ∗ is the critical cross sectional area, m2 , and subscript a and superscript * indicate nozzle outlet and critical section parameters, respectively. 2.1. Friction losses In Eq. (1), the only variable parameter is the gas dynamics function z(a ). The value of z(a ) at the supersonic region of the nozzle is an increasing function, departing from the value z( = 1) = 2.0 at the critical section. The gas dynamics function z(a ) reaches its theoretical maximum value when the dimensionless velocity reaches its maximum theoretical value max = √ (k + 1)/(k − 1) (for k = 1.2, max = 3.317). Note that the dimensionless velocity increases as the nozzle outlet cross sectional area increases, Fa . Under a constrained value of the nozzle outlet cross sectional area, Fa , the optimization of the vacuum thrust value consists of the determination of the maximum possible value of the gas dynamics function, which in the case of the nozzle with a rough internal surface depends on the nozzle opening half angle, . The theory of unidirectional compressible gases, for supersonic nozzles with rough surfaces, considers simultaneously two actions, i.e.: (i) geometric action (divergent nozzle), and (ii) friction action [15]. In order to determine the influence of friction on supersonic flow in the conic divergent nozzle, it is necessary to solve a system of differential equations that connects the velocity coefficient , with the static pressure, p (Pa), and
the nozzle cross sectional area, F, as follows [15]: 1 k 1 d () − 2 , (2) = 2 k + 1 4 tg dF¯ − 1 F¯ 2 p¯ 1 2k dp¯ = k + 1 2 − 1 F¯ () dF¯ k 2 1 × () − , k + 1 4 tg
(3)
where F¯ = F /F ∗ is the dimensionless cross sectional area; () = 1 − ((k − 1)/(k + 1))2 is a gas dynamics function; = 4f is the hydraulic resistance coefficient; f is the nozzle internal surface friction factor [16], which is taken as a nozzle parameter in the present treatment; is the conic nozzle opening half angle; p¯ = p/p ∗ is the dimensionless pressure, and p ∗ is the critical pressure, Pa. From the analysis of Eqs. (2) and (3), it follows that in the rough surface nozzle there exists a minimum opening half angle min , under which the flow does not accelerate, since d/dF¯ = 0, and the pressure does not reduce, in spite of the increase in the nozzle cross sectional area with respect to nozzle length, i.e., dp/d ¯ F¯ = 0. The numerical value of min is determined by setting the expression between brackets in Eqs. (2) and (3) to zero. Therefore tg min =
k 2 . 4 k + 1 ()
(4)
From Eq. (4), it is seen that for each combination of values of the velocity coefficient , and of the resistance coefficient there is a corresponding value of the nozzle minimum opening half angle min . At this value min , the gas adiabatic expansion work is completely absorbed by the friction forces work. Hence, the increase in gas velocity will be null, dV = 0. As a result, when the gas reaches critical parameters, an additional increase in gas velocity is possible only when the nozzle opening half angle is greater than the minimum value, i.e., > min , for some combination of and . Note also that the absolute value of min increases as the nozzle surface roughness increases. 2.2. Divergence losses Along with the roughness losses included in the gas dynamics function, non-axial divergence thrust losses have to be accounted for in the nozzle. It is a usual practice to account for the many types of nozzle thrust losses
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
869
outlet section
critical section
V flow direction F∗ α Fa Fig. 1. Axisymmetric nozzle schematic diagram.
() =
Paxial , P
1 + cos . 2
75 60
0.07 = ξ 0.05 0.03 0.01
45 30 15 0 1
3 3.317= λmax
2 λ
Fig. 2. The effect of surface roughness on minimum nozzle opening half angle with respect to dimensionless velocity.
(5) 3. Results and discussion
where () < 1 is one of possibly several correction coefficients to be applied to the ideal gross momentum thrust. The divergence losses are typically small for most nozzles (on the order of 1%), but may represent a sizeable fraction of the total nozzle performance loss and may be quite significant in nozzles with large degrees of divergence. The divergence losses of simple axisymmetric nozzles were originally analytically predicted by Malina [17]. Uniform, attached, non-interfering inviscid flow is the simplifying assumption for the analysis. Utilizing Eq. (5), vector calculus and using trigonometric relations, the final result is the familiar divergence relation for the common simple axisymmetric supersonic nozzle [17], as follows: () =
90
αmin(°)
by applying correction factors to the ideal performance characteristics of an ideal nozzle. The velocity coefficient, for example, is a factor which, when multiplied by the ideal thermodynamic velocity, accounts for some of the viscous losses in the nozzle. Similarly, correction factors for the non-axial flow of various divergent jets may be analytically derived by expressing the ratio of the momentum of all integrated axial components of the diverging flow across the nozzle exit to the momentum of the flow of an ideal nozzle where all of the exit flow is axial. Therefore, for non-ideal thrust calculations, a geometry-dependent expression for a nozzle divergence thrust loss coefficient is defined as
(6)
The values of min computed with Eq. (4) for several values of dimensionless velocity and of the hydraulic resistance coefficient (roughness), for k = 1.2, are shown in Fig. 2. From that, it is seen that the values of min increase as the flow dimensionless velocity, , and the resistance coefficient, , increase. When the dimensionless velocity reaches its maximum value, √ max = (k + 1)/(k − 1) (for k =1.2 and max =3.317), the nozzle opening half angle is the right angle, i.e., 90◦ . The determination of min has a clear practical meaning. For designing a nozzle, from given dimensionless velocity values, and resistance coefficients, the value of min determines a lower bound for a region of possible values of the nozzle opening half angles. As a result of several physical aspects, it can be concluded that for a conic supersonic nozzle, there must be an optimal nozzle opening half angle opt , under a constrained value of the nozzle outlet dimensionless
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
P˜ =
pc
Paxial ∗ F (2/k + 1)1/(k−1)
1.45
1.4
~ P
1.35 αopt 1.3
1.25 4
(7)
The objective function to be maximized is therefore the dimensionless group P˜ = z(a ) () given by Eq. (7), i.e., the product of the gas dynamics function by the thrust loss coefficient. The objective function exhibits a maximum with respect to nozzle opening half angle . Fig. 3 shows the behavior of the dimensionless group P˜ versus , with a clear maximum for all tested values of the resistance coefficient , for a fixed value of the nozzle outlet cross sectional area F¯a = 12. The analysis of Fig. 3 proceeds noting that in the conic nozzle with a smooth surface, = 0, the increase of the nozzle opening half angle causes thrust reduction, by increasing divergence losses according to Eq. (6), due to the increase in the radial velocity component at the nozzle outlet. The increase in the resistance coefficient, > 0, obviously decreases the thrust compared to the smooth nozzle. However all curves show maxima corresponding to optimal values of the nozzle opening half angle . As the nozzle internal surface roughness increases, those maxima become sharper and
8
12
16
20
24
α (°) Fig. 3. The maximization of thrust with respect to nozzle opening half angle for a fixed nozzle outlet area, F¯a = 12, and several values of nozzle internal surface roughness.
1.5
20
Fa=12 αopt ~ Pmax
1.4
10
1.3 0
= z(a )().
ξ=0 0.01 0.02 0.04 0.06
0.02
0.04
αopt (°)
cross sectional area F¯a = Fa /F ∗ , that results in maximum thrust. The existence of the thrust maximum is explained by analyzing the two extremes, i.e., small and large : (i) under constrained F¯a , the reduction of the nozzle opening half angle, , results in the increase of the area of the nozzle internal walls, which therefore causes an increase in friction losses, and (ii) on the other hand, with the increase in , friction losses decrease, but the thrust loss, (), increases due to the increase in the radial velocity component at the nozzle outlet. The optimal value of the nozzle opening half angle is determined by using the system of equations formulated by Eqs. (2) and (3). The values of were varied from 0◦ to 24◦ and the resistance coefficient assumed the values = 0, 0.01, 0.02, 0.04, and 0.06. The outlet dimensionless (relative) area assumed the values F¯a =2, 3, 6, 12, and 23. For each selected set (, , F¯a ), starting from initial conditions at the critical cross sectional area, i.e., at F¯ = 1, = 1 and p¯ = 1, the solution for and p¯ is accurately marched in the flow direction until the specified nozzle outlet area, F¯a is achieved, determining a and p¯ a . The system of equations is integrated with respect to F¯ explicitly using an adaptive time-step 4th–5th order Runge–Kutta method [18]. The F¯ -step is adjusted automatically according to the local truncation error, which is kept below a specified tolerance (10−6 in this study). A dimensionless axial thrust is defined by combining Eqs. (1) and (5) as follows:
~ Pmax
870
0 0.06
ξ Fig. 4. The behavior of the maximized thrust and optimal nozzle opening half angle with respect to nozzle internal surface roughness, for F¯a = 12.
shift to the right, towards higher values of , in order to better balance friction losses and thrust reduction. Fig. 4 shows the behavior of the maximized thrust, P˜max , and the corresponding optimal nozzle opening half angle opt with respect to the resistance coefficient, , for a fixed nozzle outlet area, F¯a =12. The maximized thrust decreases monotonically and the optimized nozzle opening half angle increases as the resistance coefficient (roughness) increases. The objective function, P˜ , is plotted in Fig. 5 with respect to the variation of the nozzle opening half angle , for several values of the dimensionless nozzle outlet cross sectional area, F¯a , and for a fixed value of the resistance coefficient, = 0.06. Again, clear maxima are obtained, for all tested F¯a values. Overall, it is seen that P˜ (or thrust) increases as the nozzle outlet cross
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
slower rate with respect to the increase in nozzle outlet area. From the point of view of practical engineering design, it is also important to observe the ‘robustness’ of the optimum found, i.e., that the variation of the optimized nozzle opening half angle is of only 2◦ for 6 F¯a < 23, therefore, for that F¯a -range, opt ∼ 15◦ for maximum thrust.
1.4 23=Fa 12 1.3 6 ~ P
αopt 3
1.2
4. Conclusions
2 1.1 4
8
12
16
20
24
α (°) Fig. 5. The maximization of thrust with respect to nozzle opening half angle for a fixed nozzle internal surface roughness, = 0.06, and several values of nozzle outlet area.
1.5
20
ξ=0.06
10
1.3 ~ Pmax
αopt (°)
~ Pmax
αopt
0
1.1 0
871
8
16
24
Fa Fig. 6. The behavior of the maximized thrust and optimal nozzle opening half angle with respect to nozzle outlet area, for = 0.06.
sectional area, F¯a , increases, for all nozzle opening half angles, . This is due to the increase in gas flow velocity. Fig. 6 shows the behavior of the maximized thrust, P˜max , and the corresponding optimal nozzle opening half angle opt with respect to the nozzle outlet area, F¯a , for a fixed value of the resistance coefficient, = 0.06. The maximized thrust and the optimized nozzle opening half angle increase monotonically as the nozzle outlet area, F¯a , increases. As the nozzle outlet area, F¯a , increases, the nozzle length increases, and the optimal nozzle opening half angle opt gradually increases as well. It is interesting to observe that the maximized thrust increases significantly with nozzle outlet area for 2 F¯a < 16, approximately. From that point on, the slope reduces and the maximized thrust increases at a
In this paper, a theoretical and numerical study was conducted to analyze hydraulic losses in conic supersonic nozzles. First the effect of the velocity and roughness on nozzle geometry (opening half angle) was investigated. Next, an optimization study was conducted to seek the possibility of finding optimal configurations (geometry) of conic supersonic nozzles for minimum thrust loss. The optimal nozzle geometry (half angle) was presented graphically, as a function of roughness (hydraulic resistance coefficient) and nozzle outlet cross sectional area. The key conclusions of engineering interest from this study are summarized as follows: 1. The existence of minimum values of the nozzle opening half angle was established for the conic nozzle, under which it is impossible for an additional increase of the gas flow velocity. That angle slowly increases as the dimensionless flow velocity and the resistance coefficient (roughness) increase. As the dimensionless flow velocity approaches its maximum value, the minimum nozzle opening half angle abruptly increases. The practical meaning of this finding is that for designing a nozzle, from given dimensionless velocity values, and resistance coefficients, the value of min determines a lower bound for a region of possible values of the nozzle opening half angles. 2. It was demonstrated that for each value of the resistance coefficient , for a fixed nozzle outlet cross sectional area, F¯a , there is an optimal nozzle opening half angle, opt , such that maximum thrust is obtained, and 3. It was also shown that for each value of the nozzle outlet cross sectional area, F¯a , for a fixed resistance coefficient , there is an optimal nozzle opening half angle, opt , such that maximum thrust is obtained. For practical engineering design, it was observed that opt ∼ 15◦ is a ‘robust’ optimum for 6 F¯a < 23, therefore, a considerably large range of variation of nozzle outlet cross sectional area.
872
D. Vlassov et al. / Acta Astronautica 61 (2007) 866 – 872
Appendix Notation cp specific heat of gas at constant pressure, J kg−1 K −1 cv specific heat of gas at constant volume, J kg−1 K −1 f friction factor F cross sectional area F¯ dimensionless cross sectional area, F (F ∗ )−1 k adiabatic ratio cp /cv p static pressure, Pa pc combustion chamber pressure, Pa p¯ dimensionless pressure, p(p ∗ )−1 P vaccum thrust, N P˜ dimensionless axial thrust, Eq. (7) R gas constant, J kg−1 K −1 Tc stagnation temperature in the combustion chamber, K V gas velocity, m s−1 z() gas dynamics function, + 1/ Greek symbols ()
nozzle opening half angle dimensionless velocity, V (V ∗ )−1 hydraulic resistance coefficient, 4f gas dynamics function, 1 − ((k − 1)/(k + 1))2 nozzledivergence thrust loss coefficient
Subscripts a axial min max
nozzle outlet axial direction minimum maximum
Superscript *
critical section (throat)
References [1] D. Vlassov, G.I. Mazing, About losses in rough surface nozzles, Journal of Machinery Construction (in Russian) 9 (1987) 81–85.
[2] J.D. Hoffman, Approximate analysis of nonisentropic flow in conical nozzles, Journal of Spacecraft and Rockets 6 (11) (1969) 1329–1337. [3] C. Hurst, A. Schulz, S. Wittig, Comparison of calculated and measured heat-transfer coefficients for transonic and supersonic boundary-layer flows, Journal of Turbomachinery —Transactions of the ASME 117 (2) (1995) 248–254. [4] M. Zaccaria, D. Ristic, B. Lakshminarayana, Three-dimensional flowfield in a turbine nozzle passage, Journal of Propulsion and Power 12 (5) (1996) 974–983. [5] A. Demetriades, Stabilization of a nozzle boundary layer by local surface heating, AIAA Journal 34 (12) (1996) 2490–2495. [6] A. Brandstein, Y. Levy, Laser propulsion system for space vehicles, Journal of Propulsion and Power 14 (2) (1998) 261–269. [7] V.V. Semyonov, Wave drag of streamlined corrugated surface of jet nozzle and flight vehicle, Izvestiya Vysshikh Uchebnykh Zavedenii Aviatsionaya Tekhnica 4 (2000) 18–22. [8] T. Ito, K. Fujii, A.K. Hayashi, Computations of axisymmetric plug-nozzle flowfields: Flow structures and thrust performance, Journal of Propulsion and Power 18 (2) (2002) 254–260. [9] C.M. Byington, J.D. Hoffman, Effects of base pressure on conical thrust nozzle optimization, Journal of Spacecraft and Rockets 7 (3) (1970) 380–386. [10] A.L. Addy, A.R. White, Optimization of drag minimums including effects of flow separation, Journal of Engineering for Industry—Transactions of the ASME 95(1), 360–364. [11] C. George, G. Candler, R. Young, E. Pfender, J. Heberlein, Nozzle optimization for dissociated species transport in low pressure plasma chemical vapor deposition, Plasma Chemistry and Plasma Processing 16 (1) (1996) S43–S56 (Suppl. S.). [12] J. Periaux, H.Q. Chen, B. Mantel, M. Sefrioui, H.T. Sui, Combining game theory and genetic algorithms with application to DDM-nozzle optimization problems, Finite Elements in Analysis and Design 37 (5) (2001) 417–429. [13] N.V. Taylor, C.M. Hempsell, Throat flow modeling of expansion deflection nozzles, JBIS—Journal of the British Interplanetary Society 57 (7–8) (2004) 242–250. [14] V.M. Kudriavtsev, Fundamental and Calculation Basis of Liquid Fuel Rocket Engines (in Russian), Superior School Publishers, Moscow, 1975, pp. 77–80. [15] A.I. Leontiev, Gas Dynamics (in Russian), MHTS Publishers, Moscow, 1997, pp. 259–279. [16] A. Bejan, Convection Heat Transfer, second ed. Wiley, New York, 1995, p. 345. [17] F. Malina, Characteristics of the rocket motor unit based on the theory of perfect gases, Journal of Franklin Institute 230 (1940) 433. [18] D. Kincaid, W. Cheney, Numerical Analysis, first ed., Wadsworth, Belmont, CA, 1991 (Chapter 8).