Investigation of Dual-Vortical-Flow Hybrid Rocket Engine without Flame Holding Mechanism
A 250 kgf thrust hybrid rocket engine was designed, tested, and verified in this work. Due to the injection and flow pattern of this engine, this engine was named dual-vortical-flow engine. This propulsion system uses N2O as oxidizer and HDPE as fuel. This engine was numerically investigated using a CFD tool that can handle reacting flow with finite-rate chemistry and coupled with the real-fluid model. The engine was further verified via a hot-fire test for 12 s. The ground Isp of the engine was 232 s and 221 s for numerical and hot-fire tests, respectively. An oscillation frequency with an order of 100 Hz was observed in both numerical and hot-fire tests with less than 5% of pressure oscillation. Swirling pattern on the fuel surface was also observed in both numerical and hot-fire test, which proves that this swirling dual-vortical-flow engine works exactly as designed. The averaged regression rate of the fuel surface was found to be 0.6~0.8 mm/s at the surface of disk walls and 1.5~1.7 mm/s at the surface of central core of the fuel grain.
Hybrid rocket propulsion has attracted tremendous attention in the past decade due to its distinct characteristics and performance as compared to liquid and/or solid rocket propulsion. The advantages of hybrid rocket propulsion include the following [1, 2]: (1) extremely high safety because of separation of oxidizer and fuel storage and solid form of fuel minimizing fuel explosion; (2) good cost-effectiveness attributing to reduced complex plumbing and valve system; (3) good throttle capability similar to liquid rocket engines; and (4) highly green and environmental friendly combustion technology which allows various choices of fuel and oxidizer. Despite the benefits mentioned above, there are some issues that a hybrid rocket engine design must take into account before it becomes a useful rocket propulsion technology. These well-known issues are mostly due to the inherent characteristics of hybrid rocket engine (HRE), which are described next.
These issues include (1) O/F ratio shift during combustion due to varying burning area and regression rate ; (2) limited total operating duration, which affects the size (diameter or length) of the hybrid rocket engine; (3) low combustion efficiency due to the nature of diffusion flame as compared to premixed flame; and (4) various choices of fuel and oxidizer [2, 4]. When designing a new type of hybrid rocket engine, one needs to take the above issues into account based on the specific mission requirements.
There are many different types of designs of HRE nowadays. For example, Nagata et al.  proposed an HRE called CAMUI (cascaded multistage impinging-jet). Fuel grain of this HRE has small ports along the engine axis, similar to the multiport design that increases the fuel burning area. The innovative part of this engine was the cascading feature of the ports. This creates large turbulence intensity as the flow impinges at the surface of the next fuel grain which boosts mixing efficiency while increasing the reaction area. A second example was proposed by Knuth et al. . This HRE came up with the idea of coaxial, corotating vortex flow field engine called Orbitec Hybrid. This HRE injects the oxidizer at the rear end of the engine in the tangential direction. After injection, the oxidizer flows to the front part of the engine along the outer surface then back to the nozzle with swirling motion. A similar engine was also studied by Wall . This engine increases the fuel regression rate by applying swirl injection , in which this design also enhances mixing by introducing bidirectional axial flow. More information on alternative designs on HRE could be found in Haag’s study . All these studies are designed to meet specific requirement missions, such as large thrust for short period of time (boosters) and small thrust for long operation (cruising). Lai et al.  recently proposed a highly efficient dual-vortical-flow (DVF) engine. The overall performance of the proposed configuration may give a very stable and high efficiency combustion and thrust. That specific configuration includes a flame holding mechanism implemented inside the high-temperature combustion chamber, which is very technically challenging or even impractical for application purposes. In this study, we would like to present a similar but easy to fabricate HRE without the flame holding mechanism for sounding rocket flight mission based on numerical and experimental investigation.
This engine consists of two counter-rotating flow reacting zones perpendicular to the engine axis with four tangential oxidizer injectors each. The fuel grain disks are connected by a central port to the nozzle, as shown in Figures 1(a) and 1(b). The main objective of this design is to provide a relatively stable thrust throughout the entire operating period. With this design, the burning surface could maintain a nearly constant combustion area throughout the operation. The counter-rotating flow maximizes the mixing and combustion efficiency with possible roll control of the rocket. Furthermore, this DVF HRE has a very small aspect ratio (L/D~1), which is favorable for gimbal-based thrust vector control (TVC) if needed.
2. Research Methods
This study initially analyzed the described engine numerically using the well-known computational fluid dynamics (CFD) technique  and was verified via hot-fire tests. The numerical tool used in this work was UNIC-UNS , which is to be described later. Figure 1(a) shows the schematic diagram of the design, in which the engine was designed to fit into a casing of 266 mm in diameter and 148 mm in length (not including the convergent-divergent nozzle). In this design, the reacting zone (blue) is surrounded by fuel grain (yellow) which regresses while combustion occurs. The arrangement of fuel grain automatically prevents the casing from the flame during operation before it is burned out. Due to the pyrolysis of the fuel grain during combustion, some of the dimensions of the engine change over time, namely, Pd, Dd, and Gap in Figure 1(a). The dimensions for the numerical investigation are summarized in column “Cases A1~A3” of Table 1. The pyrolysis rate of the fuel grain is a function of many physical properties such as temperature, flow pattern, turbulence intensity, and oxidizer mass flux, to name a few.
2.1. Numerical Analysis
The engine model was simplified and solved numerically by computational fluid dynamics (CFD) method according to the governing equations, namely, continuity equation (mass conservation equation), species conservation equations, Navier-Stokes equations (momentum conservation equation), and energy conservation equation with various well-known physical models such as HBMS (Hirschfelder, Beuler, McGee, and Sutton) real-fluid model  and extended k-ε turbulence model . These equations and models were discretized using the cell-centered finite-volume method, parallelized using MPI protocol and had been applied successfully to similar problems [10, 12]. For the reacting species, the finite rate reaction model (a.k.a. Arrhenius reaction model) was implemented. Together with a simplified set of species and reaction path , the UNIC-UNS code handles the species conservation equation. In general, the simulations were performed using a time step of 5 × 10−6 s with boundary conditions summarized in Table 2. These computations are performed using 64 cores of processors on the IBM 1350 PC cluster at the National Center for High-Performance Computing (NCHC) of Taiwan, in which a single node consists of 4 cores with 3.0 GHz of CPU for each core and 16 GB of RAM for each node.
For the CFD model, we performed a series of grid convergence tests using 5.03, 2.34, and 1.63 million cells which were labeled as Case A1, A2, and A3, respectively. We have compared three simulated quasi-steady-state data which include thrust, mass flow rate (ṁ), and Isp. Figure 2 illustrates the corresponding convergence history, in which all the cases were calculated long enough to reach a quasi-steady-state where the abovementioned values are nearly constant. The thrust oscillation caused by the pressure oscillation will be explained later in Results and Discussion. For all three cases, the simulated mass flow rates are essentially the same as 1.05 kg/s. But the thrust (or Isp) of Case A3 (coarse mesh) is slightly higher than the A1 (fine mesh) and A2 (medium mesh) which are ~254 kgf (or 234 s) and ~246 kgf (or 232 s), respectively. With these results, one can summarize that the resolution of Case A2 is good enough for engine design purpose, further analysis, and discussion. The corresponding hot-fire tests are described next.
2.2. Hot-Fire Test Setup and Analysis
Figure 3 shows the schematic sketch of the DVF engine design in detail. For hot-fire test purpose, the chamber body, external plumbing, and both bulkheads were made of 304 stainless steel. The nozzle was fabricated using high-density antioxidation graphite. The insulators are made of silicone rubber, EPDM, or ceramic depending on the location. The fuel grains were manufactured using off-the-shelf high-density polyethylene (HDPE)  with a density of 0.945 g/cm3. Photographs of all manufactured components are shown in Figure 4(a) and are assembled into a DVF engine as shown in Figure 4(b).
Figure 5 shows the test setup of this hot-fire test. For this system, the N2 cylinder was filled to exceed 120 atm and was connected to the top of the specially prepared oxidizer (N2O) running tank with a regulator set as 57 atm. A pressure transducer was also mounted to the top of the running tank to monitor the tank pressure during operation. The bottom opening of the running tank was connected to each of the injectors of the engine with a plumbing system including a main valve followed by a flow distributor and some pressure transducers. The DVF engine was then mounted on a horizontal thrust stand. The pipe size used before the flow distribution system was 3/4 inch (~19 mm) and was split into eight 1/4-inch (~6.35 mm) tube before injection into the combustion chamber. The valve used in this system was a traditional ball valve with nominal diameter 19 mm controlled by a pneumatic valve.
A reliable electrical ground support equipment (EGSE) is also required to accomplish the task. To acquire the experimental data, we mounted a series of sensors on the engine system: these sensors were managed on a data acquisition system (DAQ) based on National Instruments Corporation’s (NI) products. The sensors used in this test were pressure transducers and load cells. These sensors output 0~10 VDC according to the physical quantity measured. The pressure transducers were JPT-131 series provided by Jetek Electronics Co. Ltd. . The pressure transducer can measure 0~100 bar gauge pressure with the accuracy of +/−0.5%. For the thrust of the engine and mass difference of the supply tank, S-type load cell with the desired load range provided by Sensolink  was used. These sensors were wired to the cRIO 9074 of NI using ordinary signal cables. Then the cRIO 9074 is connected to a PC using Ethernet cable, which is suitable for long-distance monitoring and controlling as shown in Figure 6. The programming was done using the software called LabVIEW, which is also provided by NI. This platform makes the programming process simple and easy. The obtained data is stored in the PC at the rate of 1 kHz. These data are then postprocessed and analyzed after the hot-fire test.
3. Results and Discussion
3.1. Simulation Results
We have performed the calculation of thermodynamic equilibrium reactions for N2O  (oxidizer) and C2H4 (fuel) using NASA CEA online [17–19]. The chamber pressure was set to 38 atm with an O/F ratio of 4.2, and the area expansion ratio of the nozzle was 2.56. For this case, the resulting outlet pressure 2.85 atm was underexpanded for sea level operation. The optimal (theoretical) Isp, C, and Cf value for this case was 236 s, 1759 m/s, and 1.32, respectively. Figure 7(a) shows the summary of the simulated equilibrium mass fractions of all species. The major species include N2, CO, H2O, CO2, and H2 with the mass fraction of 0.5136, 0.3572, 0.0664, 0.0417, and 0.0200, respectively. These five species sum up to 0.9989 of the composition. Figure 7(b) shows the simulated mass fractions by UNIC-UNS CFD code. The major species are identical to those of CEA, but the mass fractions are slightly different. The mass fraction of N2, CO, H2O, CO2, and H2 is 0.4715, 0.2924, 0.0713, 0.0554, and 0.0121, respectively. In addition to the above species, there are other species that sum up to 0.0973, which are mostly radical species of combustion. The main reason for this difference was that for CEA the species were calculated assuming thermodynamic equilibrium while, for UNIC, they were obtained using the finite-rate chemistry. The majority of the mixture is N2, which comes from the direct decomposition of N2O oxidizer under high combustion temperature exceeding 1640°C. Due to the low O/F ratio which leads to incomplete combustion, the second dominant species is CO. The formation and distribution of the species will be further discussed later in this paper.
Figure 8(a) shows the instantaneous sliced pressure distribution in the engine in the quasi-steady-state after ignition of combustion. The pressure inside the engine ranges from 30 to 38 atm, in which the maximum pressure is located near to the injectors (circumference of the disks) and the lowest pressure is distributed near the nozzle as expected. Noticeably, the pressure in the major central core (between disk 1 and disk 2) is low with a value in the range 32-33 atm due to strong swirling motion. This strong swirling central core region disappears as the counter flows meet somewhere downstream of disk 2. As the flow continues to flow towards the nozzle, the flow accelerates and the pressure drops quickly. In addition, Figure 8(b) shows the same instantaneous sliced distribution of Mach number in the engine. Gas flow inside the engine is subsonic, which is accelerated to sonic speed at the throat of the nozzle and is further accelerated in the diverging part.
The performance of HRE relies heavily on the combustion efficiency, which temperature of combustion can be considered as a good indicator. Figure 9(a) shows the corresponding instantaneous sliced temperature distribution. Near the injection regions where N2O are injected, the temperature is at room temperature (300 K), which serves as a natural cooling mechanism for preventing the injectors from melting by high combustion temperature. As the N2O stream is injected into the engine, the high combustion temperature causes N2O to decompose directly. This decomposition reaction is an exothermic reaction that helps to sustain the combustion. The main products of N2O decomposition are N2 and O2 with some related species in radical form. Distributions of some critical species such as O2 and OH radical are shown in Figures 9(b) and 9(c), respectively. Therefore, we can observe that in Figure 9(b), massive O2 are formed and quickly disappear just before the high-temperature region. In hydrocarbon combustion, flame location or highly reacting region generally consists of abundant OH radical generally. This shows that the abundant OH radical distribution corresponds to the high-temperature region very well.
Figure 10 shows the distributions of mass fraction of the five other major species (N2, CO, H2O, CO2, and H2) in the engine. In Figure 10(a), the mass fraction of N2 is about 0.6 in the bulk region of the two major combustion disk chambers due to the direct decomposition of N2O. It decreases rapidly in the very thin region near the fuel grain surface because of the dilution from fuel vapor due to pyrolytic reaction. As the flow goes further downstream, the mass fraction of N2 decreases to about 0.45 before entering the divergent part of the nozzle. Figures 10(b), 10(c), and 10(d) show the distribution of mass fraction of H2O, CO2, and CO, respectively, which represents the major products of complete combustion of hydrocarbon compound with oxygen. Note that the distribution of CO2 in the disk chambers (Figure 10(c)) is approximately in consistent with that of OH (Figure 9(c)) which corresponds to the flame position where strong chemical reaction occurs. The results show the combustion changes from fuel-lean, stoichiometric, and fuel-rich reactions along the radial direction from injection to central core region. The concentrations of H2O and CO2 increase along the radial direction inwards because of reactions between decomposed O2 and fuel grain vapor and reach maximal values (H2O: 0.09, CO2: 0.15) somewhere before the flowing into the central core. In addition, they are brought to concentrate near the fuel grain surface of the central core because of the centrifugal force of swirling flow. Due to the lack of “O” species (Figure 9(b)) (more fuel-rich), both the concentrations of CO2 and H2O decrease after entering the central port. With the same reason, a large amount of syngas (i.e., CO and H2) is formed because of the high temperature before (in the disk chamber) and near the wall in the central port region (Figure 9(a)). These species (CO and H2) can be further reacted in this central port region as the flow moves further downstream, especially the case of H2 (Figure 10(e)).
We have found that the thrust (or pressure, not shown) oscillates as a function of time in Figure 2. A detailed analysis on this oscillating phenomena was performed in this study. Figure 11 illustrates five instantaneous temperature distribution in the middle sliced sections of both disk combustion chambers at five different times of simulation. The “spikes” of the contour rotate clockwise and counterclockwise at disk chambers 1 and 2, respectively. The oscillation frequency based on numerical simulations is about 200 Hz which should be due to the lack of flame holding mechanisms in the reacting region near the disks. This oscillation is more distinct as compared to our previous work of DVF HRE with flame holders . The reason without using flame holders as explained earlier is that the implementation of the flame holding mechanism becomes very difficult or even impossible due to the high temperature (exceeds 3000 K) in the reacting region with highly vibrating environment and the fuel shape (internal geometry) changes as it operates. In spite of the lack of flame holders, the performance of the engine is considered to be fairly stable with the maximal oscillation amplitude of 2 kg over 245 kg, which is less than 1%.
Figure 12 shows the stream traces of the flow patterns at the surface of the fuel grain (0.1 mm above the surface) of Case A2 (Figure 12). The injected flow streams revolve in the disk combustion chambers around the central port before entering the central port. This provides a relatively longer flow path for the flow (relative to straight radial injection) which greatly increases the residence time for the combustion reactions to take place. For a specific case, the stream trace marked with the blue arrow in Figure 12 was being observed. We define the stream trace entering the swirl pattern when it reaches the radial position where R is at 95% of Rmax (about 5 mm from the circumference, indicated by the long red arrow) and exit when R is at 25% of Rmax (about 5 mm larger than the central port, indicated by the short red arrow). This figure shows that the flow revolved about 180 degrees starting from injection until it entered the central port. The path that the flow takes in this situation is about twice the length of the radius of the disk. As shown in Figure 12, the gap between disk chamber walls is small and the averaged flow speed is larger. This also indicates that the tangential momentum is larger; therefore, the flow revolves almost 180 degrees before reaching the central port. As the gap grows wider, the cross-sectional area increases and the tangential momentum becomes smaller. This will be further compared and discussed with the hot-fire test results.
3.2. Hot-Fire Test Results
After the rocket engine is set up following Figure 5 with all connections carefully checked, we have followed the following procedure to perform the hot-fire tests. A snapshot of the hot-fire test during combustion is shown in Figure 13, in which the exhausting plume is slightly underexpanded because of clear further expansion leaving the lip of the nozzle. Figure 14 shows the measured thrust and pressures at many locations for a typical run. A pyrograin in the engine is used to ignite the engine at t = −4 s. After the pyrograin burns out, main valve opens at t = 0 s to allow the N2O oxidizer to flow into the chamber. A delay of 0.5 s could be observed based on the rise of measured thrust and pressure data, which is caused by the speed of valve opening. The liquid N2O flow depletes at t = 9.5 s, and finally, the valve totally closes at t = 12 s to shut down the engine, during which both thrust and pressures decrease rapidly. Note that the N2O becomes a gaseous state and flows into chamber between t = 9.5 and t = 12 s. The thrust is relatively stable in the range of 240~245 kgf in the period of 1.25~9.5 s. The thrust decreases almost exponentially from t = 9.5 s to the end. The running tank pressure starts at 57 atm and decays slowly to about 50 atm at t = 9.5 s due to cooling effect caused by thermal expansion of N2O flow during operation. Both disk chamber pressures show clear oscillations (even different) probably due to the lack of flame holding mechanism in the engine which coincides with the simulation results. However, the oscillation frequency is only 4 Hz based on the measured pressure data in Figure 14, which is much lower than the experimentally observed 100 Hz of oscillation using 600 fps high-speed camera most probably caused by the limitation of the setup of pressure sensor that may damp out the high-frequency component. In addition, the simulated oscillation frequency is ~200 Hz which is roughly consistent with the measured ~100 Hz. The difference, however, requires further investigation.
The swirl pattern is also observed on the surface of the fuel grain after hot-fire test. Figure 15 shows the photo of cavity wall of chamber disk 1 after hot-fire test. The flow pattern was indicated by the red arrows. Similar to the simulated case, the two red arrows in Figure 15 indicate the entering and exiting of the flow into the chamber disk. The revolved angle of this case is about 80 degrees which is a lot smaller than the simulated Case A2. The main reason was that the gap after firing (~26 mm) is a lot larger than the initial case (10 mm). As the gap increases, the ratio of tangential momentum and radial momentum decreases. Therefore, the angle revolved by the flow is expected to decrease before entering the central port.
The fuel surface contour was measured by a bridge-type 3D Coordinate Measuring Machine (CMM) (Model PIONEER, Hexagon Manufacturing Intelligence). We have scanned three fuel grain surfaces using an automatic mode. Figure 16 shows the averaged regression rates at different radial positions by taking average from 6 to 12 scans per radius along different radial directions, considering minimal measured positions (0, 45, and 48 mm for fuel grain 1, 2, and 3, resp.) to the outer radial position of 100 mm. The averaged regression rate at the disk surface is 0.6 to 0.8 mm/s. The regression rate at the center region of grain 1 (central port) is two times the value of those at the disk surface. This highly regressed region may be attributed to the long light of sight of radiation from the combustion flame to the exhausting plume along the axis. In addition, very high regression rates are also observed at the wall of the central port which is about 1.5 to 1.7 mm/s, which should be caused by the very strong swirling in this region that promotes the pyrolysis of the fuel grain.
3.3. Comparison of Simulation with Test Data
Table 3 summarizes the numerically simulated cases with hot-fire test data of the 250 kgf class DVF HRE. For the numerical results, there are two specific cases being discussed. The NASA CEA case assumed 0-D chemical equilibrium condition which we can be considered as the theoretical results, and the UNIC-UNS simulated Case A2 using CFD finite-rate model. For CEA test case, the initial state of the engine is calculated. The resulting O/F ratio, ground Isp, and C are 4.2, 236 s, and 1759 m/s, respectively. The corresponding Cf is calculated to be 1.32 using the nozzle area expansion ratio of 2.56.
For the CFD simulated cases, the simulated Case A2 represents the initial state of this work. The simulated Case A2 has a 98% of Isp efficiency (232 s) and more than 100% of the C value of 1814 m/s which exceeds that (1759 m/s) of CEA. For this work, the location Pc used is somewhere near the injector where the pressure sensor of the experimental model could be mounted during tests. Due to the use of this pressure, the C value was higher than the one from CEA. As a result, the Cf value of Case A2 is not as high as the one from CEA. This can be easily observed by .
For the hot-fire test results, the pressure and thrust are obtained as a function of time but the mass difference (especially the fuel part) could only be measured after the test due to the limitation of the facilities. Therefore, the mass flow rate of the hot-fire test was only available in terms of time averaged values. Due to the reasons above, this work could only obtain the averaged value of Isp and C of the hot-fire test. But to obtain the value of Cf, the instantaneous value can be calculated properly by Cf = Finstant/(Pinstant × Athroat). The averaged Isp was calculated by integrating the measured force and divide it by total mass difference. The C value was then calculated with Pc taken as time averaged value. As for the Cf, it was calculated using the instantaneous force and chamber pressure. The averaged Isp obtained by the hot-fire test in this work was 221, about 93.6% of the theoretical value. The averaged C and Cf are 1565 m/s and 1.39, respectively. Since the value of Cf was available as a function of time, we observed an increase of Cf from 1.21 to 1.51, which is mainly due to the decrease of chamber pressure as the test proceeds. As the chamber pressure decreases, the underexpanded flow shifts towards optimal criteria, and the Cf increases. Though the main objective of this DVF HRE was to provide a constant and stable thrust for a specific flight mission, it was surprising, that to us, that the thrust remains almost the same throughout the test period. This was probably due to O/F ratio shift, mixing efficiency, and a lot more reasons that compensate with each other. However, this definitely requires further investigation in the near future.
This work proposes a DVF HRE without flame holding mechanism for possible sounding rocket application. The length-to-diameter ratio of the engine is only about 1, which is different from conventional lengthy HREs. This design was simulated considering geometrical configurations at the initial state. The resulting simulated Isp at initial state is 232 s with very high combustion efficiency as compared to that calculated by NASA CEA (236 s). A hot-fire test based on proposed design was performed for 12 s by measuring the thrust and pressures at many locations. The maximal thrust of 245 kgf was measured to be relatively constant with a value of ~240 kgf and an averaged ground Isp of 221 s. Measured regression rates are in the range of 0.6–0.8 mm/s at the walls of the two disk chambers and 1.5–1.7 mm/s at the walls of the central port region due to strong vortex motion of the hot gases. In addition, the central end wall, located furthest from the nozzle, also possesses very high regression rate, probably due to high thermal radiation and highly turbulent flow field. The O/F ratio of this specific test was 4.2 which is relatively low compared with the value for optimal Isp from various references . An oscillation frequency of ~100 Hz was observed in numerical simulation, while ~200 Hz in hot-fire test, which definitely requires further investigation. These observed or measured instabilities may be caused by the fact that no physical flame holders were used in the chamber. Despite with the presence of these instabilities, this DVF engine design still shows fairly stable thrust, which should satisfy the use in sounding rocket application.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by the Advanced Rocket Research Center (ARRC), National Chiao Tung University of Taiwan, and Ministry of Science and Technology of Taiwan (MOST) through Grant nos. 102-2627-E-009-001, 103-2221-E-009-136-MY2, and 105-2221-E-009-066-MY2 and technical assistance of Taiwan Innovative Space Inc. of Taiwan.
G. P. Sutton and O. Biblarz, Rocket Propulsion Elements, JOHN WILEY & SONS, INC., 7th edition, 2001.
M. J. Chiaverini and K. K. Kuo, Fundamentals of Hybrid Rocket Combustion and Propulsion, AIAA, 2006.
W. J. Knuth, D. J. Gramer, M. J. Chiaverini, J. A. Sauer, R. H. Whitesands, and R. A. Dill, “Preliminary CFD analysis of the vortex hybrid rocket chamber and nozzle flow field,” in 34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Cleveland, OH, USA, July 13–15, 1998.View at: Publisher Site | Google Scholar
N. Wall, Characterisation of Multiple Concentric Vortices in Hybrid Rocket Combustion Chambers, Ph.D. Thesis, Department of Mechanical Engineering, University of Sheffield, South Yorkshire, England, 2013.
G. S. Haag, Alternative Geometry Hybrid Rockets for Spacecraft Orbit Transfer, Ph.D. Dissertation, School of Electronic Engineering, Information Technology and Mathematics, University of Surrey, Surrey, United Kingdom, 2001.
J. Anderson, E. Dick, G. Degrez et al., Computational Fluid Dynamics an Introduction, Springer-Verlag, Berlin Heidelberg, 2009.
S. Gordon and B. J. McBride, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications, I Analysis, NASA RP-1311, 1994.
B. J. McBride and S. Gordon, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications, II Users Manual and Program Description, NASA RP-1311, 1996.