NUMERICAL STUDY OF THE AERODYNAMIC CHARACTERISTICS OF
A WAVE DISC ENGINE
By
Guangwei Sun

A DISSERTATION
Submitted to
Michigan State University
In partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2012

ABSTRACT
NUMERICAL STUDY OF THE AERODYNAMIC CHARACTERISTICS OF
A WAVE DISC ENGINE
By
Guangwei Sun
The demand for mobile power generation and low-cost durable and yet efficient and low
emission engines is ever increasing. This work focuses on the aerodynamic characteristics of a
novel combustion engine – the wave disc engine. This new core technology uses a turbo
combustion shock wave technique to efficiently convert fuel sources to mechanical power and
subsequently electrical power. In such engines, sudden closing of the flow passage is used to
generate shock wave that can enhance combustion and allow for compression and expansion of
the operating gases. This new engine concept combines the advantages of (a) high efficiency
confined combustion comparable to automotive internal combustion engines (IC-E) and
cutting-edge pulse detonation engines (PDE), with (b) the high power density and low
maintenance of continuous flow gas turbines (GT), but (c) at a much lower unit cost due to its
physical simplicity and compactness.
This work performs a detailed thermodynamic analysis of a Wave Disc Engine by taking a
single channel as the control volume as well as a full rotor. The analysis gives the general
expression of the net work done by the flow and the thermal efficiency in terms of the
mass-average properties, and states that under some certain conditions the operation can achieve
the efficiency of Humphrey cycle. Two one-dimensional numerical codes were created using

TVD-MacCormack scheme and Gottlieb-Turkel scheme to provide an estimation of the ports
timing for multi-dimensional design. A simple arc was chosen as the channel shape of the
baseline engine. After some geometric adjustments, the configuration of the baseline design was
determined. However, its efficiency is very low. Energy losses were investigated by a 2D
simulation of a single channel. The results showed that the significant outgoing kinetic energy
and the under-expansion of the exhaust gas were the main sources of energy loss. To improve the
performance of the energy extraction by the curved channels, the mechanism of the torque
generation by the expansion wave was studied and the principle for the channel shape design was
obtained. Several channel shapes were proposed and tested. Eventually the convergent channel
was chosen as the best design and implemented on a full engine. The results showed that the
efficiency of the engine increased from 2.2% to 3.3% by using the convergent channel. Then the
effects of the peak pressure in the channel and of the tip speed on the engine power and
efficiency was investigated. Two methods to reuse the energy from the exhaust gas were
proposed – adding a return channel or an external turbine. It was proved by the simulations that
both methods can improve the efficiency and adding an external turbine is more effective. At last,
various plans to prevent the leakage going into the inlets were compared and the jagged wall
design was found to be the best.

Copyright by
GUANGWEI SUN
2012

TABLE OF CONTENTS
LIST OF TABLES ....................................................................................................................... VII
LIST OF FIGURES ...................................................................................................................... IX
NOMENCLATURE ................................................................................................................... XIV
CHAPTER 1. INTRODUCTION................................................................................................. 1
1.1. Wave Rotor Technology ................................................................................................... 1
1.1.1.
Wave Rotor ............................................................................................................... 1
1.1.2.
Wave Engine ............................................................................................................. 3
1.1.3.
Wave Combustor ....................................................................................................... 4
1.2. Working Principle of the WDE ........................................................................................ 5
1.3. Thermodynamic Cycle Analysis ...................................................................................... 6
1.3.1.
Single Channel Analysis ........................................................................................... 7
1.3.2.
Full Rotor Analysis ................................................................................................. 13
CHAPTER 2. 1D SIMULATION............................................................................................... 16
2.1. Numerical Scheme ......................................................................................................... 17
2.1.1.
Gottlieb-Turkel Scheme .......................................................................................... 18
2.1.2.
TVD-MacCormack Scheme.................................................................................... 21
2.1.3.
Boundary Conditions .............................................................................................. 24
2.2. Determination of the Ports Timing ................................................................................. 25
2.3. Test Case and Results ..................................................................................................... 26
2.1. Particle Tracking ............................................................................................................ 29
CHAPTER 3. 2D SIMULATION............................................................................................... 30
3.1. Optimization Parameters ................................................................................................ 30
3.2. Simulation methodology ................................................................................................ 31
3.3. Optimization for a Given Rotor ..................................................................................... 33
3.4. Baseline Design .............................................................................................................. 34
CHAPTER 4. METHODS TO IMPROVE THE EFFICIENCY ............................................... 37
4.1. Investigation into the Energy Usage .............................................................................. 37
4.1.1.
1D Investigation ...................................................................................................... 37
4.1.2.
2D Investigation ...................................................................................................... 37
4.2. Investigation into the Channel Shape ............................................................................. 42
4.2.1.
Torque Generation Mechanism for Unsteady Flow ................................................ 43
4.2.2.
Effect of the Outlet Opening ................................................................................... 49
4.3. Implementation of the New Shape ................................................................................. 55
V

4.4. Effect of the Peak Pressure............................................................................................. 56
4.4.1.
1D Investigation ...................................................................................................... 56
4.4.2.
2D Investigation ...................................................................................................... 57
4.4.3.
Full Engine Test ...................................................................................................... 57
4.5. The Effect of the Blade Tip Speed ................................................................................. 58
4.5.1.
Single Channel Investigation .................................................................................. 58
4.5.2.
Full Engine Investigation ........................................................................................ 62
4.6. Reuse of the Energy from the Exhaust Gas.................................................................... 63
4.6.1.
Return Channel ....................................................................................................... 63
4.6.2.
Two-Stage Rotor ..................................................................................................... 64
4.6.3.
Configuration of the Guide Vanes and the Turbine Blades ..................................... 65
4.6.4.
Scalability ............................................................................................................... 66
CHAPTER 5. METHODS TO PREVENT BACKFIRE............................................................ 68
5.1. Simulation Methodology ................................................................................................ 68
5.2. Various Backfire Prevention Designs ............................................................................. 73
5.2.1.
No prevention design .............................................................................................. 73
5.2.2.
Plan A: U channel ................................................................................................... 73
5.2.3.
Plan B: modified U channel .................................................................................... 74
5.2.4.
Plan D: leakage pipe ............................................................................................... 76
5.2.5.
Plan E: jagged housing wall .................................................................................... 76
CHAPTER 6. CONCLUSIONS AND FUTURE WORK .......................................................... 79
6.1. Summery and Conclusions ............................................................................................. 79
6.2. Suggestions for Future Work .......................................................................................... 81
APPENDIX….. ............................................................................................................................. 84
REFERENCES. .......................................................................................................................... 157

VI

LIST OF TABLES
Table 2-1: Dimension of the computational domain and the working condition of the test case for
the validation of the 1D code ........................................................................................................ 27
Table 2-2: Ports timing obtained by GT scheme and TVDM scheme .......................................... 28
Table 3-1: Symbols in Figure 3-1 and their meaning ................................................................... 31
Table 3-2: Working conditions for the baseline design ................................................................. 32
Table 3-3: Final geometric parameters of the baseline WDE ....................................................... 35
Table 3-4: Power and efficiency of the baseline design................................................................ 36
Table 4-1: Results of the 2D simulation for the analysis of the energy loss ................................. 40
Table 4-2: Work generated in each process in the sixth cycle ...................................................... 44
Table 4-3: Results of the unsteady simulations for C1 and C2 channel ....................................... 45
Table 4-4: General principle of torque generation by travelling expansion wave ........................ 47
Table 4-5: Results of the single channel simulations for C3, J1 and J2 channels......................... 49
Table 4-6: Outlet openings, expansion durations and efficiencies of the single channel
simulations for J2, C3, MJ2 and MC3 channels ........................................................................... 50
Table 4-7: Blade angles of the convergent channel ...................................................................... 50
Table 4-8: Results of the single channel simulations for the convergent channel ........................ 51
Table 4-9: Results of the simulation with patching method for the baseline design .................... 56
Table 4-10: Results of the WDE with the convergent channels.................................................... 56
Table 4-11: Results of the initial pressure tests ............................................................................. 57
Table 4-12: Results of the simulation with 2100 K patching temperature for the WDE with C3
channels......................................................................................................................................... 58
Table 4-13: Results of the simulation with 1200 K patching temperature for the WDE with C3
channels......................................................................................................................................... 58

VII

Table 4-14: Results of the rotor dimension tests ........................................................................... 60
Table 4-15: Results of the rotational speed tests ........................................................................... 61
Table 4-16: Comparison of the ports timing of different rotational speeds .................................. 62
Table 4-17: Results of the WDE with 30 cm OD ......................................................................... 63
Table 4-18: Results of the WDE with a return channel ................................................................ 63
Table 4-19: Results of the two-stage engine ................................................................................. 64
Table 4-20: Results of the two-stage engine with a traditional external turbine........................... 66
Table 4-21: Results of the scaled engine....................................................................................... 67
Table 5-1: Working conditions of the simulations in CHAPTER 5 .............................................. 70
Table 5-2: cp at different temperatures at 1 atm [19] .................................................................... 71
Table 5-3: Results of various backfire prevention designs ........................................................... 78

VIII

LIST OF FIGURES
Figure 1: Schematic model of a wave rotor inside the pressure wave supercharger (PWS) with its
development stages (soldered, cast with single of channels, cast with double row of channels)
and actual NASA four port wave rotor drum (right). For interpretation of the references to color
in this and all other figures, the reader is referred to the electronic version of this dissertation. . 84
Figure 2: Porting and flow scheme for a through-flow (left) and reverse flow wave rotor (taken
from [55]) ...................................................................................................................................... 85
Figure 3: Assembly of successfully working wave engine with wave rotor (left) and rear and
front stator plates [3] ..................................................................................................................... 85
Figure 4: ABB test rig (left) and wave rotor combustor cross section (right) [4] ......................... 86
Figure 5: By using shockwaves for internal work transmission, the Wave Disk Engine virtually
eliminates the hardware (shown in gray) in gas turbines and internal combustion engines that
serve to transmit expansion work internally to perform the needed inlet gas compression. ........ 86
Figure 6: Schematic of an internal combustion WDE .................................................................. 87
Figure 7: T-s diagram of the mass-average flow in a WDE .......................................................... 87
Figure 8: Wave diagram obtained from the TVDM scheme ......................................................... 88
Figure 9: Profiles of pressure, density, velocity and temperature for TVDM scheme (black) and
GT scheme (blue) at 5.44e-5 s and its corresponding location in the wave diagram ................... 89
Figure 10: The profiles of pressure, density, velocity and temperature from GT scheme (blue) and
TVDM scheme (black) at 7e-5 s and its corresponding location in the wave diagram ................ 90
Figure 11: Profiles of pressure, density, velocity and temperature from (b) GT scheme (blue) and
TVDM scheme (black) at 2.3e-4 s and its corresponding location in the wave diagram (a)........ 91
Figure 12: Profiles of pressure, density, velocity and temperature from (b) GT scheme (blue) and
TVDM scheme (black) at 5.33e-4 s and its corresponding location in the wave diagram (a)...... 92
Figure 13: T-s diagram (upper) and p-v diagram (lower) for one fluid particle ........................... 93
Figure 14: Geometric parameters of a WDE ................................................................................ 94
Figure 15: Heater for ignition ....................................................................................................... 95

IX

Figure 16: Backflow into the air inlet ........................................................................................... 95
Figure 17: Relative velocity vector plot (a) and contour of mass fraction of methane (b) ........... 96
Figure 18: Absolute velocity vector plot near the inlet (a) and outlet (b) ..................................... 96
Figure 19: Backfire into the fuel inlet ........................................................................................... 97
Figure 20: Incomplete combustion ............................................................................................... 98
Figure 21: (a) temperature contour, (b) mass fraction of methane contour and (c) relative velocity
vector plot for the baseline design ................................................................................................ 99
Figure 22: Power history in the second cycle ............................................................................. 100
Figure 23: T-s diagrams of 10 particles in the channel ............................................................... 101
Figure 24: Geometry for the 2D investigation ............................................................................ 102
Figure 25: Histories of the outlet enthalpy (left) and outlet pressure (right) of the 2D simulation
..................................................................................................................................................... 102
Figure 26: Geometry of the C2 channel ...................................................................................... 103
Figure 27: Pressure contours (Pa) in the curved channels .......................................................... 103
Figure 28: History of the moment coefficient on the channel walls ........................................... 104
Figure 29: Geometry of the baseline shape................................................................................. 104
Figure 30: Profiles of a shock wave (top) and a planar acoustic wave (bottom) through a
two-dimensional right-angle bend, numbers 1 - 6 are representative wave front configurations
within the bend (Edwards et al [12], Kentfield [23]) .................................................................. 105
Figure 31: Sequence of the pressure contours when the channel just opens .............................. 106
Figure 32: Sequence of the pressure contours during the incident expansion wave propagates 107
Figure 33: Sequence of the pressure contours when the reflecting wave interferes with the
incident wave .............................................................................................................................. 107
Figure 34: Pressure contours as the expansion wave propagates in a straight channel .............. 108
Figure 35: Pressure contour of the C2 channel ........................................................................... 108
Figure 36: Geometry of the C3 channel ...................................................................................... 108
X

Figure 37: Geometries of J1 (left) and J2 (right) channels ......................................................... 109
Figure 38: Efficiency vs. expansion duration for C1, C2, C3, J1 and J2 channels ..................... 109
Figure 39: Outlet opening (circle radius) .................................................................................... 109
Figure 40: Expansion duration vs. outlet opening .......................................................................110
Figure 41: Geometries of the MJ2(left) and MC3(right) channels ..............................................110
Figure 42: Geometry of the convergent channel .......................................................................... 111
Figure 43: Torque histories of five short expansion cases (upper) and three long expansion cases
(lower) ..........................................................................................................................................112
Figure 44: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the propagation of the incident
expansion wave ............................................................................................................................113
Figure 45: Pressure contour (top) and pressure profiles (bottom) near the trailing edge of the
channel walls when the channel outlet is not fully open .............................................................115
Figure 46: Pressure contour (top), pressure profiles (middle) and Mach number contour near the
trailing edge when the channel outlet is fully open .....................................................................116
Figure 47: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the propagation of the reflected
expansion wave ............................................................................................................................118
Figure 48: Relative velocity vectors showing the back flow near the inlet end ..........................119
Figure 49: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the positive pressure gradient caused
by the back flow near the inlet .................................................................................................... 120
Figure 50: Pressure contour in the beginning of the simulation ................................................. 121
Figure 51: Torque histories of the first three cycles.................................................................... 122
Figure 52: Geometry of the single rotor engine with convergent channels ................................ 123
Figure 53: Histories of torque of the first two cycles ................................................................. 124
Figure 54: Profile of pressure in case 1 when the incident expansion propagates ..................... 125

XI

Figure 55: Profile of pressure in case 2 when the incident expansion propagates ..................... 125
Figure 56: WDE with C3 channels ............................................................................................. 126
Figure 57: Velocity triangle at the trailing edge of a channel wall ............................................. 126
Figure 58: Efficiency vs. outer radius ......................................................................................... 127
Figure 59: Efficiency vs. RPM ................................................................................................... 127
Figure 60: Geometry of the engine with 30 cm OD ................................................................... 128
Figure 61: WDE with a return channel ....................................................................................... 129
Figure 62: Configuration of a two-stage engine ......................................................................... 130
Figure 63: Pressure contour in the external turbine .................................................................... 131
Figure 64: Streamlines in the guide vane passage ...................................................................... 132
Figure 65: Configuration of a second generation engine ............................................................ 133
Figure 66: Pressure contour in the blade passages of the second generation engine .................. 134
Figure 67: Streamlines in the guide vane passage of the second generation engine .................. 135
Figure 68: Temperature contour of the no prevention design at 2e-5 s ...................................... 136
Figure 69: Thermodynamic process of the hot leakage gas ........................................................ 137
Figure 70: Energy patching area (red area) ................................................................................. 138
Figure 71: Variation of cp with T of the experimental data and the fitted curve ........................ 139
Figure 72: Temperature contour at the end of the third cycle from the pressure patching method
for the no prevention design ....................................................................................................... 140
Figure 73: Idea of the U channel................................................................................................. 141
Figure 74: Temperature contour at the end of the third cycle from the pressure patching method
for the U channel design ............................................................................................................. 142
Figure 75: Velocity vector plot sequence of the left U channel during one channel period ....... 142
Figure 76: Velocity vector plot sequence of the right U channel during one channel period ..... 144

XII

Figure 77: Geometry of the MU channel .................................................................................... 146
Figure 78: Velocity vector plot sequence of the left MU channel during one channel period.... 146
Figure 79: Velocity vector plot sequence of the right MU channel during one channel period . 148
Figure 80: Temperature contour at the end of the third cycle from the pressure patching method
for the MU channel design.......................................................................................................... 150
Figure 81: Comparison between the pressure in the MU channels (top) and U channels (bottom)
..................................................................................................................................................... 151
Figure 82: Configuration of a WDE with the leakage pipes ....................................................... 152
Figure 83: Temperature contour at the end of the third cycle from the pressure patching method
for the leakage pipe design ......................................................................................................... 153
Figure 84: Geometry of the jagged housing wall........................................................................ 154
Figure 85: Temperature contour at the end of the third cycle from the pressure patching method
for the jagged wall design ........................................................................................................... 155
Figure 86: Using EGR to generate a shock wave and ignition ................................................... 156

XIII

NOMENCLATURE
A

area

c

absolute velocity

cv

constant volume specific heat

cp

constant pressure specific heat

Cr

courant number

E

total energy in the control volume

h

specific enthalpy

H

enthalpy

k

kinetic energy

𝑚̇

mass flow rate

mch

mass in one channel

N

node number

Nch

channel number of one cycle

p

pressure

Q

heat release

s

specific entropy

T

temperature

Tcycle

cycle period

Tch

channel period
XIV

u

specific internal energy

W

work done by the fluid

α

ports wall angle

β

blade angle

ρ

density

η

thermal efficiency

γ

specific heat ratio

ω

angular velocity

τ

torque on the channel walls

Δ

discrepancy

Subscripts and Superscripts
1

left boundary node

2

the beginning of the constant volume combustion

3

the end of the constant volume combustion

CV

control volume

c

compression process

e

expansion process

in

inflow

mean

mean value

net

net value
XV

new

new time step

old

old time step

origin

origin point of characteristic line

out

outflow

r

radial component

ref

reference value

s

scavenging process

tot

total property

�
”
θ

̇

mass-average property
time rate of property
per unit area
tangential component

XVI

CHAPTER 1. INTRODUCTION
1.1. Wave Rotor Technology
The idea of the energy exchange between two mediums without any separation goes back to the
beginning of the 20th century. In its second decade, along the longitudinal axis perforated drum,
a channeled rotor was patented by the German engineer Burghard - a machine delivering an
uninterrupted mass flow of the pressurized air [3]. As the unsteady flow theory, a necessity for
the development of a usable machine, was not developed until the 1920’s and 1930’s, Burghard’s
invention did not succeed in becoming an available device until 1940’s when the Brown Boveri
(BBC, today ABB) turbocharger engineer Seippel designed a pressure exchanger, also called the
wave rotor, as an air compressor of a gas turbine used as the propulsion of an experimental
locomotive [3].

1.1.1.

Wave Rotor

The wave rotor is a non-steady flow device that uses shock waves to pressurize fluids by
transferring energy from a high-pressure flow to a low-pressure flow in narrow channels. For a
gas turbine engine application, the wave rotor employs the hot, high pressure exhaust gas from
the combustion chamber to generate a shock wave that compresses the cooler, lower pressure air
received from the compressor. This results in an increase in both temperature and pressure of the
air entering the combustion chamber, allowing for a higher overall pressure ratio for the entire
cycle for the fixed turbine inlet temperature. Such pressure exchange wave rotor effectively

1

combines the steady-state traditional turbomachinery with unsteady, compressible gas flow
principles, to allow for higher cycle efficiency.
The wave rotor, shown in Figure 1, consists of a series of axial channels in a rotating drum,
driven either by a motor or by V belt direct from the internal combustion (IC) engine crank shaft.
The wave rotor was firstly successfully commercialized as an IC engine supercharger by Mazda
in 1980’s and 90’s in serial production of diesel passenger cars [52]. The air inlet represents the
low-pressure air received from the compressor (in case of the gas turbine) or from the ambient
air (in case of the IC engine supercharger), exhaust inlet is the high-pressure gas (exhaust) from
the combustor or from cylinders of the IC engine, exhaust outlet is the low-pressure gas
delivered to the turbine from the wave rotor or to the low pressure part of the IC engine exhaust
manifold, and air inlet is the high pressure-air delivered to the combustor or to the cylinders of
the IC engine from the wave rotor. The rotor channels are placed between two end plates with
control ports. The number of ports can vary. Usually, as two cycles per rotor revolution are
realized, enabling equal thermal stress, four ports at each of the plates are placed. Each port is
designed to expose the channels to the working fluids at a specific shaft angle and for a specific
duration. Shock and expansion waves are initiated inside the channels by pressure differences,
caused by the port opening and closing. Because the channels are exposed to both hot and cool
gases, wave rotors can be naturally self-cooled. Due to the pre-expansion of the burned gases in
the wave rotor besides the rotor self-cooling feature, the combustor can operate at higher
temperatures without raising the turbine inlet temperature. This is especially important in
applications where the temperature is limited by material constraints. Schematic diagrams of the
2

flow through the channels both for the through flow and reverse flow configurations are shown
in Figure 2.

1.1.2.

Wave Engine

The wave engine consists of a single wave rotor, used for compression of the inlet air and
expansion of the exhaust gas, while producing power. It includes a single steady flow gas turbine
type combustor.
In the UK of the mid-1950s, the Ruston-Hornsby Turbine Company, manufacturer of
diesel engines and industrial gas turbines, supported the construction and testing of a different
kind of wave rotor designed by Pearson [3]. This unique wave rotor, known as the wave engine,
has helical axially arranged channels that change the direction of the gas flows producing shaft
work similar to a conventional turbine blade. The rotor has a 230 mm diameter and a 76 mm
length. The engine worked successfully for several hundred hours in a wide range of operating
conditions (from 3000 to 18,000 rpm) without variable porting, and produced up to 26 kW at its
design point with a cycle peak temperature of 1070 K and a thermal efficiency of around 10%.
Better performance seemed possible with improvements to overcome leakage and incomplete
scavenging. Even though large number of seizures occurred due to unexpectedly large axial
displacement of the rotor, no engine damage resulted. The design of the engine was based on
wave diagrams using the method of characteristics (Figure 2).The sealing and bearings were
carefully adapted, considering rotor thermal expansion. Eventually, the engine was wrecked due
to over-speeding from an improperly connected fuel line, and the project was canceled.

3

1.1.3.

Wave Combustor

The wave combustor is a wave rotor with the internal combustion which can produce a pressure
rise, in comparison to the steady flow combustion which is always accompanied by a pressure
loss.
ABB in collaboration with the Swiss Federal Institute of Technology Zurich commenced
the design of a full wave rotor with 36 combustor axial channels in 1992 [4].

Each channel had

165 mm length and 15×15 mm cross section (Figure 4). The 200 mm inner-diameter rotor was
driven by an electric motor capable of up to 5000 rpm. The spark-plug electric ignition was used
for the start up, and the auto injection by the hot gas injection was used for the continuous
operation. To start up, two spark plugs with sparking frequency of 3000 Hz located at each end
of the rotor and for each cycle would ignite the combustible mixture in the channels from both
sides simultaneously. Self-sustaining ignition was accomplished by employing axial jet injection
of post combustion gas from a neighboring channel. Although the fuel injection and ignition
finally worked fine, tests revealed following shortcomings 1) inhomogeneous mixture in the cell,
resulting in a slow diffusion flame, 2) maximum pressure of 9 bar, relative to expected 10–15 bar,
presumably due to leakage, which also caused premature ignition and misfiring at higher
chamber pressures, 3) excess thermal stresses on the ignition ring, 4) inadequate cantilever
single-bearing rotor support, 5) complicated and sensitive electromechanical device for
controlling leakage gap. Major remedies recommended making the system better included 1)
lead away duct for leakage gas removal, 2) rotor cooling by air, 3) two-sided rotor support, and 4)
mechanical control for thermal expansion. The prototype wave combustor operated successfully
4

until the project was concluded in 1994, due to market concerns.

1.2. Working Principle of the WDE
The Wave Disc Engine (WDE) utilizes a combustion engine cycle consisting of compression,
combustion, expansion with work extraction, and heat rejection to the ambience in narrow, radial
arranged and curved channels.
The compression work is typically provided through a work generated during expansion. A
considerable part of the machinery hardware in the GT and IC-E is devoted to extracting and
transmitting the work needed for compression (gray in Figure 5). In a WDE this is mainly
performed via shock and pressure waves that are formed in the working fluid and triggered by
sudden opening and closing of the channels in the rotating disk to the inlet and outlet ports.
Because of this, the hardware of the WDE reduces to about one 1/4 to 1/3 that is necessary to
extract the brake horsepower (useful external work) shown in red in Figure 5.
Also, using shockwaves that move with sonic speed reduces the inertia of the hardware
and ensures rapid response. With the expansion at sonic speed immediately after combustion, the
residence time at high temperature is extremely short, even shorter than in a gas turbine. This
results in ultra-low NOx emissions. Furthermore, the relevant heat transfer time and areas are
extremely small – dramatically reducing heat losses.
Figure 6 shows a schematic of a WDE. The cycle begins with the sudden closing of the
outlet port. A hammer shock wave is generated by the deceleration of the exiting flow to zero
velocity and propagates towards the inlet, which in turn compresses the fresh air-fuel mixture
5

behind the shock wave. The inlet port is still open to allow for more loading. Once the shock
wave arrives at the inlet end, the inlet port closes and the mixture in the channel is ignited. After
the ignition, the constant volume combustion takes place within the channel, producing a
pressure and temperature rise during the combustion process. After the combustion is completed,
the outlet end of the channel opens to the ambience. This sudden opening of the channel creates
an expansion wave propagating towards the inlet. Torque generation is produced by the fluid
tangential momentum at the outlet (jet propulsion). Once the pressure near the inlet end of the
channel is lowered down by the expansion wave below the inlet pressure, the inlet port opens
and the fresh air-fuel mixture at the inlet pressure is drawn into the channel and flushes the
exhaust gas away. Centrifugal force on the flow helps the flushing and loading processes. When
the channel is filled up with the mixture, the outlet port closes suddenly. Then the cycle starts
over.

1.3. Thermodynamic Cycle Analysis
Researchers have performed the thermodynamic analysis of various unsteady flow devices, such
as PDE, wave rotor/combustor topped gas turbine, using either Brayton cycle or Humphrey cycle.
Weber's engine operated on the Brayton cycle (Weber [67]). Heiser and Pratt used Humphrey
cycle for PDE (Heiser and Pratt [65]). Wintenberger and Shepherd used FJ cycle to model PDEs
(Wintenberger and Shepherd [66]). Wilson and Paxson used Brayton cycle to model the wave
rotor topped gas turbine (Wilson and Paxson [67]). Kentfield modeled PDEs as a Lenoir cycle
(Kentfield [23]).
6

These analyses were based on the energy equation for steady flow or closed system. Since
the flow in the WDE is unsteady as in a PDE, not all working fluid particles undergo identical
thermodynamic processes. Besides, a WDE cannot be treated as a closed system. Thus, the
energy equation for steady flow or closed system does not apply on the WDE. A new analytical
way is presented here by integrating the following energy equation over the control volume and
the time:

dECV
=
Q + m in htot ,in − m out htot ,out − W
dt

(1.1)

where all physical quantities vary with time due to the unsteady nature of the flow, and W
indicates the work done by the flow through the control volume. If the value of W is negative, the
flow consumes energy; if the value is positive, the energy is extracted from the flow. In this
analysis, the leakage is neglected and all walls are assumed to be adiabatic.

1.3.1.

Single Channel Analysis

When the WDE operates in a periodic mode, each channel experiences the same working process.
Thus, any channel can represent all channels. In this analysis, a single channel was taken as the
control volume.
From the closing of the outlet (t1) to the closing of the inlet (t2), the channel is filled up with
fresh fuel-air mixture which is compressed by the hammer shock wave during this process.
Integrating Eq. (1.1) over this duration, we obtain
t2

t2
dECV
=
∫ dt dt ∫ m in htot ,in dt − Wc
t1

t1

7

(1.2)

Rearranging Eq. (1.2), we obtain the compression work done by the fluid in the channel in this
compression process, which is
t2

t2

t1

t1

=
Wc ∫ m in htot ,in dt −

=

∫

dECV
dt
dt

t2

∫ m in htot ,in dt + E1 − E2

(1.3)

t1

From the closing of the inlet (t2) to the opening of the outlet (t3), the constant-volume
combustion takes place in the channel. Integrating Eq. (1.1) from t2 to t3, we obtain the heat
addition to the channel, which is
Q = E3 − E2 = 
∫∫∫ ρ3u3dV − ∫∫∫ ρ2u2dV
CV

(1.4)

CV

From the opening of the outlet (t3) to the opening of the inlet (t4), the high pressure gas
discharges to the ambience and the expansion wave propagates towards the inlet. Integrating Eq.
(1.1) from t3 to t4, we have
t4

∫

t3

t

4
dECV
dt =
− ∫ m out htot ,out dt − We
dt

(1.5)

t3

Rearranging Eq. (1.5) we obtain the work done by the gas in the channel, which is
t4

t4

t3

t3

We =
− ∫ m out htot ,out dt −

∫

dECV
dt
dt

t4

=
− ∫ m out htot ,out dt + E3 − E4
t3

8

(1.6)

From the opening of the inlet (t4) to the closing of the outlet (t5/t1), the fresh mixture is
sucked into the channel and the residual exhaust gas is scavenged. At the end of this process, the
channel returns to the beginning point of the working cycle. Integrating Eq. (1.1) from t4 to t5,
we have
t5

∫

t4

t5
t5
dECV
dt =∫ m in htot ,in dt − ∫ m out htot ,out dt − Ws
dt
t4

(1.7)

t4

Rearranging Eq. (1.7) we obtain the work done by the flow in this process:
t5

t5

t5

t4

t4

t4

t5

t5

Ws =∫ m in htot ,in dt − ∫ m out htot ,out dt −

=

∫

dECV
dt
dt

∫ m in htot ,in dt − ∫ m out htot ,out dt + E4 − E5

t4

t4

where E1 = E5 due to the periodicity.
Adding all three expressions of the work done by the fluid, we obtain
Wnet = Wc + We + Ws


 t2
  t4

= ∫ m in htot ,in dt + E1 − E2  +  − ∫ m out htot ,out dt + E3 − E4 


 
 t1
  t3

t5
 t5


+ ∫ m in htot ,in dt − ∫ m out htot ,out dt + E4 − E5 


t4
 t4


=

t2

t4

t5

t5

t1

t3

t4

t4

∫ m in htot ,in dt + E3 − E2 − ∫ m out htot ,out dt + ∫ m in htot ,in dt − ∫ m out htot ,out dt

9

(1.8)

t2

t5

t5

t1

t4

t3

t2

t5

t5

=
Q + ∫ m in htot ,in dt + ∫ m in htot ,in dt − ∫ m out htot ,out dt

=
Q + ∫ m in htot ,in dt + ∫ m in htot ,in dt − ∫ m out htot ,out dt
t1

t4

(1.9)

t3

Due to the periodicity,
t2

=
∫ m in htot ,in dt
t1

t2 +Tcycle

=
∫ m in htot ,in dt

t2 +Tcycle

∫

t1 +Tcycle

m in htot ,in dt

(1.10)

t5

Substituting Eq. (1.4), Eq. (1.10) into Eq. (1.9), we have
t2 +Tcycle

t5

t5

t5

t4

t3

Wnet =
Q+
∫ m in htot ,in dt + ∫ m in htot ,in dt − ∫ m out htot ,out dt

=
∫∫∫ ρ3u3dV − 
∫∫∫ ρ2u2dV +
CV

t2 +Tcycle

∫

m in htot ,in dt − ∫ m out htot ,out dt (1.11)

t4

CV

t5

t3

Eq. (1.11) gives the general form of the work output of the channel during the cycle.
The general thermal efficiency is

η
=

t +T
t
Q + ∫ 2 cycle m in htot ,in dt − ∫ 5m out htot ,out dt

Wnet
=
Q

t4

t3

Q
t2 +Tcycle

t5

=

m in htot ,in dt
∫ t3m out htot ,out dt − ∫ t4
1−
∫∫∫ ρ3u3dV − 
∫∫∫ ρ2u2dV

CV

(1.12)

CV

During the combustion, the mass in the channel consists of two parts: 1) the existing mass

10

in the channel at t1 which is

∫∫∫ ρ1dV ; and 2) the new mass entering the channel from t1 to t2

CV

t2

which is

∫ m in dt .

t1

At t1 the existing mass in the channel is equal to the mass that enters from t4 to t5 due to the
periodicity, i.e.,
t5

∫∫∫ ρ1dV = ∫ m in dt

CV

(1.13)

t4

Since the exhaust gas in the channels is fully scavenged by the fresh mixture and no mixture
overflows out of the channels when the outlet is closed, the mass contained in the closed channel
is incoming mass during the scavenging and compression processes. Combined with the mass
conservation, we have
t2

t5

t2

t1

t4

t1

mch = 
∫∫∫ ρ2dV = 
∫∫∫ ρ3dV = 
∫∫∫ ρ1dV + ∫ m in dt = ∫ m in dt + ∫ m in dt
CV

CV

CV

t2 +Tcycle

=

∫

m in=
dt m=
in

t4

t5

dt
∫ m out=

mout

(1.14)

t3

Defining mass-average specific internal energy and static temperature as

u2
=

∫∫∫ ρ2u2dV


CV
=

∫∫∫ ρ2dV

CV

11

cv T 2

(1.15)

∫∫∫ ρ3u3dV


u3
=

CV
=

∫∫∫ ρ3dV


(1.16)

cv T 3

CV

Defining mass-average total specific enthalpy and total temperature as
t2 +Tcycle

∫ t4

m in htot ,in dt
= c p T tot ,in
t2 +Tcycle
m in dt
∫

=
htot ,in

(1.17)

t4

t2 +Tcycle

=
htot ,out

∫ t4

m out htot ,out dt
= c p T tot ,out
t2 +Tcycle
m out dt
∫

(1.18)

t4

Substituting Eq. (1.15) - (1.18) into Eq. (1.12), we have

h

(

c p T tot ,out − T tot ,in

−h
u3 − u 2

tot , out
tot ,in
=
η=
1−
1−

(

cv T 3 − T 2

)

)=
T tot ,out − T tot ,in
1− γ
T3 −T 2

(1.19)

Eq. (1.19) gives the thermal efficiency in terms of mass-average temperature.
Further assumptions can be made. If the incoming kinetic energy is equal to the outgoing
kinetic energy, i.e.,
t2 +Tcycle

∫

t4

t5
2
2
cin
cout


min
dt = ∫ mout
dt
2
2

(1.20)

t3

Eq. (1.12) becomes
t2 +Tcycle

t5

η=

m in hin dt
∫ t3m out hout dt − ∫ t4
1−
∫∫∫ ρ3u3dV − 
∫∫∫ ρ2u2dV

CV

(1.21)

CV

Then all total properties can be replaced with the corresponding static properties, and we
have

12

η = 1− γ

1.3.2.

T out − T in
T3 −T 2

(1.22)

Full Rotor Analysis

Now take the entire rotor as the control volume. When the WDE operates in a periodic mode, the
period is the time that it takes for a channel rotates to the position of the next channel. To
distinguish with the cycle period, this period is named the channel period Tch. Applying Eq. (1.1)
to the rotor, we have

∫
Tch

dECV
dt =
dt

∫

 +
Qdt

Tch

∫ ∫∫

∫ ∫∫

′′ htot ,in dAdt −
m in

Tch inlet

′′ htot ,out dAdt −
m out

Tch outlet

∫

 =0
Wdt

(1.23)

Tch

Where ECV , 𝑄̇ and 𝑊̇ varies with time, and the total enthalpy flow rate not only varies with
time, but also has a distribution along the inlet and outlet surface.
Rearranging Eq. (1.23), we have the work output:

∫
Tch

 +
 =
Wdt
∫ Qdt
∫ ∫∫ m in′′ htot ,in dAdt −
Tch

Tch inlet

∫ ∫∫

′′ htot ,out dAdt
m out

(1.24)

Tch outlet

Since the exhaust gas in the channels is fully scavenged by the fresh mixture and no mixture
overflows out of the channels when the outlet is closed, the mass in the rotor is the mass coming
into the rotor during one cycle period. Thus, the mass in any closed channel is the incoming mass
during one channel period, and we have

mch =

∫ ∫∫

′′ dAdt
m in

(1.25)

Tch inlet

During one cycle period, the heat is released in Nch channels, meaning that in average the
heat is released in one channel during one channel period. The heat release rate is

13

) (

(

 = m u3 − u 2 = u3 − u 2
Qdt
ch

∫
Tch

) ∫ ∫∫ m in′′ dAdt

(1.26)

Tch inlet

Combining Eq. (1.24) and Eq. (1.26), we have the thermal efficiency

Wdt

∫
=
η

 +
Qdt

∫

Tch
T
= ch

∫ ∫∫

′′ htot ,in dAdt −
m in

Tch inlet


Qdt

∫

∫ ∫∫

∫

Tch


Qdt

′′ htot ,out dAdt
m out

Tch outlet

Tch

∫ ∫∫

′′ htot ,in dAdt −
m in

T inlet
= 1 + ch

∫ ∫∫

′′ htot ,out dAdt
m out

Tch outlet

(1.27)

( u3 − u 2 ) ∫ ∫∫ m in′′ dAdt
Tch inlet

Mass conservation gives

∫ ∫∫

∫ ∫∫

′′ dAdt =
m in

Tch inlet

m o′′ut dAdt

(1.28)

Tch outlet

Due to the periodicity, the total enthalpy coming through the inlet during one channel period
is equal to that coming into one channel during the scavenging and compression processes,
namely
t2 +Tcycle

∫ ∫∫

′′ htot ,in dAdt =
m in

Tch inlet

∫

m in htot ,in dt

(1.29)

t4

Combining Eq. (1.14) and Eq. (1.25), we have
t2 +Tcycle

=
mch

=
∫ m in dt
t4

∫ ∫∫
Tch inlet

Thus, Eq. (1.17) can be rewritten as

14

′′ dAdt
m in

(1.30)

∫ ∫∫
=
htot ,in

′′ htot ,in dAdt
m in

Tch inlet
= c p T tot ,in
′′ dAdt
m in
Tch inlet

(1.31)

∫ ∫∫

In the same manner Eq. (1.18) can be rewritten as

∫ ∫∫
=
htot ,out

′′ htot ,out dAdt
m out

Tch outlet
=

∫ ∫∫

′′ t dAdt
m ou

(1.32)

c p T tot ,out

Tch outlet

Substituting Eq. (1.31) and Eq. (1.32) into Eq. (1.27), we have the same form of the thermal
efficiency as Eq. (1.19)

h

−h
u3 − u 2

(

c p T tot ,out − T tot ,in

tot , out
tot ,in
=
η=
1−
1−

(

cv T 3 − T 2

)

)=
T tot ,out − T tot ,in
1− γ
T3 −T 2

(1.19)

which shows a consistency between the full rotor analysis and the single channel analysis.
In the sense of mass-average flow, all fluid particles can be considered to be of the
mass-average properties at any instance throughout their passing through the channel. The
thermodynamic cycle can be depicted on a T-s diagram (Figure 7). When they enter the rotor,
they are at the state of “in”. After being compressed by the hammer shock wave, they are at the
state of “2”. During the combustion, their state changes from “2” to “3” following the
constant-volume curve. Then the burned gas expands to the state of “out”. If the “out” state has
the identical pressure of the “in” state, the cycle that the particles undergo is a Humphrey cycle
(or Atkinson cycle).

15

CHAPTER 2. 1D SIMULATION
A key step in the design process is to determine the ports timing, i.e., when the ports should be
open or closed. When the outlet is open to the ambience, the expansion wave propagates into the
channel towards the inlet end. When reaching the inner housing wall, the expansion wave
reflects. Since there is no analytical solution for the reflected waves, a numerical method should
be used to explore the wave timing. Once the reflected expansion wave lowers the pressure at the
inlet end of the channel down to the ambient pressure, the inlet port should be open. Then the
fresh mixture is drawn into the channel. Once the channel is filled up with the mixture and
exhaust gas is cleared away from the channel, the outlet port should be closed. Then a hammer
shock wave is generated by the closing of the outlet and propagates towards the inlet. Once the
shock wave arrives at the inlet, the inlet port should be closed. Thus, the determination of the
ports timing depends on the wave pattern.
Spalding at Imperial College seems to be the first investigator using a 1D numerical method
to investigate unsteady flows inside wave rotors in early 1970s. He has formulated a numerical
procedure considering the effects of heat transfer and friction. The program was later developed
by Jonsson (Jonsson et al [23]), and it was successfully applied to pressure exchangers
calculations (Matthews [27], Azim [7], Azoury and Hai [8]). In the same time, at ABB in
Switzerland, Zehnder [69]-[71] developed and used a finite difference code for evaluation of the
flow in pressure wave machines. By initiating a wave rotor research at NASA Lewis (now Glenn)
Research Center in the late 1980s, Paxson [35][36] has developed a one-dimensional gas

16

dynamic model to calculate design geometry and off-design wave rotor performance. The model
solves the unsteady flow field in an axial passage for time-varying inlet and outlet port
conditions, using simplified models to account for losses due to gradual passage opening and
closing, viscous and heat transfer effects, leakage, and non-uniform port flow field mixing.
Recent improvement and validations have completed it as a preliminary and general design tool
Paxson [37]-[40],[41],[42]).
Fatsis and Ribaud from the French National Aerospace Research Establishment (ONERA)
have also developed a one-dimensional numerical code based on an approximate Rieman solver
taking into account viscous, thermal, and leakage losses (Fatsis and Ribaud [43]-[16]). The code
has been also compared with experimental data and applied to three-port, throughflow, and
reverse-flow configurations.
Group guided by Reizes in Australia tried to make some attempts to wave rotor technology
[59] performing several calculations of shock wave motion in rotating shock tubes. Have been
undertaken also attempts of using analytical methods for design of wave rotor cycles (Resler,
Moscari and Nalim [56]).
In 2007, Piechna created a 1D code using the Lax-Wendroff scheme taking into account the
changing cross sectional area (Piechna [52]).

2.1. Numerical Scheme
First, a 1D code was made to simulate the flow in a rotating straight channel for the preliminary
investigation. It is assumed that the flow is along the length the channel (x-direction), the
17

variables do not change along the width (y-direction), the walls are adiabatic and there is no
leakage between the channel and the housing walls. All kinds of gas in the channel were treated
as ideal gas and the viscous forces were neglected.
First, the Gottlieb-Turkel scheme was applied to solve the governing equation

∂U
∂F
=
−
+J
∂t
∂x
ρ
 
where=
U =
ρu 
 ρe 
 

 ρ 


 ρ u  is the solution vector, F=
ρc T 
 v 

(2.1)

 ρu
 2
ρu +
 ρ eu




p =



 ρu 
 2

 ρ u + p  is the flux
 ρ c Tu 
 v 



 0 


vector, J =  ρω 2 x  is the source term. In this chapter, u represents the velocity and e the
 ∂u 
− p 
 ∂x 
internal energy.
The element equations are the continuity equation, the x-momentum equation and the
energy equation.
Along with Gottlieb-Turkel scheme, fourth order artificial dissipation terms were also
explicitly added to reduce the numerical oscillations near the sudden discontinuities (Anderson
[1]). For comparison, a TVD-MacCormack (TVDM) scheme was used to solve the problem as
well. Both schemes are based on the MacCormack scheme.

2.1.1.

Gottlieb-Turkel Scheme

In 1969, MacCormack introduced a two-stage numerical scheme for compressible flows with a

18

predictor stage followed by a corrector stage (MacCormack [62]). The scheme is second-order
accurate in both space and time and is widely applicable, in part, because of its simplicity and
robustness. Gottlieb and Turkel improved the MacCormack scheme to fourth-order in space and
second-order in time. This scheme has been popular among researchers involved with highly
resolved flow fields (Gottlieb and Turkel [63]).
The Gottlieb-Turkel (GT) scheme consists of two steps, predictor step and corrector step.
The procedure of the GT scheme is described as follows.
1) Set the initial values of ρ, p, T, u and calculate the initial U and F at all grid points.
 ∂U
2) In the predictor step, the predicted time derivative 
 ∂t

t
 
 
i 

*

at the internal points is

calculated based on the known values at time level t by replacing the spatial derivative on
the right-hand of Eq. (2.1) with a fourth order 3-point backward differencing, as follows.
 ∂U

 ∂t

t *

 7 F t − 8F t + F t 

i
i −1
i − 2 =

=
−
+ J it , i 2,3, , N − 1




∆
6
x
i 


t +∆t

( )i

3) The predicted solution vector at time level t + ∆t , U
t +∆t

(U )i

=

U it

 ∂U
+ 
 ∂t

(2.2)

is calculated using Eq. (2.3).

*

t
 
t + Sit , i 2,3, , N − 1
  ∆=
i 

(2.3)

where Sit is the artificial viscosity term, calculated using Eq. (2.4).

=
Sit

Cx

Pit+1 − 2 Pit + Pit−1
Pit+1 + 2 Pit

+ Pit−1

(Uit+1 − 2Uit + Uit−1 ) ,
19

0.01 ≤ C x ≤ 0.3

(2.4)

4) Use

t +∆t

(U )i

t +∆t

to calculate the predicted values of ρ i
t +∆t

p = ρ RT to calculate pi
t +∆t

Ti

t +∆t

t +∆t

, ui

and pi

t +∆t

, and then calculate F i

t +∆t

, Ti

t +∆t

and J i

t +∆t

, ui

and use
t +∆t

based on ρ i

,

.

 ∂U
5) In the corrector step, the time derivative of the predicted solution vector, 
 ∂t
t +∆t

calculated based on the F i

t +∆t

and J i

t +∆t



i

, is

by replacing the spatial derivative on the

right-hand of Eq. (1.1) with fourth order 3-point forward differencing, as follows.

 ∂U

 ∂t

t +∆t



i

 −7 F t +∆t + 8 F t +∆t − F t +∆t
i
i +1
i+2
=
−
6∆x




t +∆t
 +=
, i 2,3, , N − 1
Ji



(2.5)

 ∂U
6) The average value of the time derivative is obtained from the arithmetic mean of 
 ∂t

t
 
 
i 

*

t +∆t

 ∂U 
and 

 ∂t i

.

 ∂U t +∆t 
1   ∂U
 =

 

2   ∂t
 ∂t i

avg


t +∆t 
t

   ∂U  =
 ， i 2,3, , N − 1

  +
i   ∂t i


*

(2.6)

7) The corrected solution vector at time level t + ∆t , U it +∆t , is calculated using Eq. (2.7).
U it +∆t=

U it

 ∂U
+ 
 ∂t

t +∆t 



i

t +∆t

， i 2,3, , N − 1

∆t + S=
i

avg

(2.7)

where
t +∆t

t +∆t

=
Si
Cx

t +∆t

Pi +1 − 2 Pi

t +∆t

+ Pi −1

 U t +∆t − 2U t +∆t + U t +∆t  ,
i
i −1 
t +∆t
t +∆t
t +∆t  i +1

Pi +1 + 2 Pi
+ Pi −1

0.01 ≤ C x ≤ 0.3 (2.8)

8) Use U it +∆t to calculate the corrected values of ρit +∆t , Tit +∆t , uit +∆t and use p = ρ RT to

20

calculate pit +∆t , and then calculate Fit +∆t and J it +∆t based on ρit +∆t , Tit +∆t , uit +∆t
and pit +∆t .
9) Go back to step (3) to run a new iteration until the final time.
In this code, the time step size is determined by the Courant-Friedrichs-Lowry (CFL)
criterion and given by

( )

t
∆t=
min=
∆tit , i 2,3, , N − 1

(2.9)

where
∆tit =
Cr

(

∆x
uit + ait

)

(2.10)

with Cr being the Courant number and set to be 0.8 in this code.

2.1.2.

TVD-MacCormack Scheme

Total-variation-diminishing (TVD) algorithm is used to eliminate the numerical oscillations in
the vicinity of discontinuities. In this section, a TVDM scheme (Dongfang Liang et al [64]) is
used. This scheme includes the addition of a five-point symmetric total variation diminishing
(TVD) term to the corrector step of the standard MacCormack scheme.
The procedure of the TVDM scheme is described as follows.
1) Set the initial values of ρ, p, T, u and calculate the initial U and F at all grid points.
 ∂U
2) In the predictor step, the predicted time derivative 
 ∂t

t
 
 
i 

*

at the internal points is

calculated based on the known values at time level t by replacing the spatial derivative on
the right-hand of Eq. (2.1) with a second order backward difference, as follows.
21

 ∂U

 ∂t

*

t
 Ft − Ft 
 
i −1 =
− i
+ J it ， i 2,3, , N − 1
  =


∆x
i 


t +∆t

( )i

3) The predicted solution vector at time level t + ∆t , U
t +∆t

(U )i
t +∆t

( )i

4) Use U

t +∆t

*

t +∆t

t +∆t

t +∆t

, and then calculate F i

, is calculated using Eq. (2.12).

t
 
 ∆t , i 2,3, , N − 1
=
i 

to calculate the predicted values of ρ i

to calculate pi
and pi

 ∂U
U it + 
=
 ∂t

(2.11)

t +∆t

, Ti

t +∆t

and J i

(2.12)

t +∆t

, ui

and use p = ρ RT

t +∆t

based on ρ i

t +∆t

, Ti

t +∆t

, ui

.

 ∂U
5) In the corrector step, the time derivative of the predicted solution vector, 
 ∂t
t +∆t

calculated based on the F i

t +∆t

and J i

t +∆t



i

, is

by replacing the spatial derivative on the

right-hand of Eq. (2.1) with forward difference, as follows.
 ∂U

 ∂t

t +∆t



i

 F t +∆t − F t +∆t
i +1
i
=
−
∆x




t +∆t
 + J=
， i 2,3, , N − 1
i



(2.13)

 ∂U
6) The average value of the time derivative is obtained from the arithmetic mean of 
 ∂t

t
 
 
i 

*

t +∆t

 ∂U 
and 

 ∂t i

.
t +∆t
t *
 ∂U t +∆t 
1   ∂U    ∂U  
=
 
 ， i 2,3, , N − 1
=
 
  +
 ∂t i  avg 2   ∂t i   ∂t i 

(2.14)

7) The corrected solution vector at time level t + ∆t , U it +∆t , is calculated using Eq. (13).

22

U it +∆t=

U it

 ∂U
+ 
 ∂t

t +∆t 



i

=

∆t + TVD
i ， i 2,3, , N − 1

avg

(2.15)

where

( ) ( )(

) ( ) ( )(
U t − U it−1 )(U it+1 − U it )
+ ( i
r =
(Uit+1 − Uit )(Uit+1 − Uit )

+
− 
+
− 
t
t
t
t


=
TVD
i G ri + G ri +1  U i +1 − U i − G ri −1 + G ri  U i − U i −1

U it − U it−1 )(U it+1 − U it )
(
r− =
(Uit − Uit−1 )(Uit − Uit−1 )

)

(2.16)

(2.17)

(2.18)

The function G() is defined as

where the flux limiter function is given as

and the variable C is

𝐺(𝑥) = 0.5𝐶[1 − 𝜑(𝑥)]

(2.19)

𝜑(𝑥) = max�0,min(2𝑥, 1)�

(2.20)

𝐶𝑟(1 − 𝐶𝑟), 𝐶𝑟 ≤ 0.5
0.25,
𝐶𝑟 > 0.5

(2.21)

𝐶=�

8) Use U it +∆t to calculate the corrected values of ρit +∆t , Tit +∆t , uit +∆t and use p = ρ RT to
calculate pit +∆t , and then calculate Fit +∆t and J it +∆t based on ρit +∆t , Tit +∆t , uit +∆t
and pit +∆t .
9) Go back to step (3) to run a new iteration until the final time.
The time step size is determined in the same way as in the previous section.

23

2.1.3.

Boundary Conditions

Three types of boundary conditions (B.C.s) were used in this code:
Inflow
The stagnation pressure and temperature are imposed; the velocity and density are
determined by Piechna's method (Piechna [52]) which is a combination of energy equation
2
ain
a12 u12
=
+
γ −1 γ −1 2

(2.22)

and the compatibility relation

a1new =
aorigin +

γ −1
2

(u1new − uorigin )

(2.23)

The static pressure is obtained from the isentropic relation
γ −1
T1new  p1new  γ
=


Ttot ,in

 ptot ,in 



(2.24)

Outflow
If the outflow is subsonic, the static pressure is imposed; temperature and velocity are
determined by a simple linear extrapolation method using the information of the interior nodes; if
the outflow is supersonic, the static pressure, temperature and velocity are all determined by the
linear extrapolation (Anderson [1]):

− p N − 2 ,?supersonic
2 p
p N =  N −1
pout ,?subsonic


(2.25)

=
u N 2u N −1 − u N − 2

(2.26)

=
TN 2TN −1 − TN − 2

(2.27)

24

Solid Wall
In the case of right closed end wall the simplest but efficient B.C.s are following:
p N = p N −1

(2.28)

u N = −u N −1

(2.29)

TN = TN −1

(2.30)

They represent the existence of the rigid wall between node N and node N-1. In the similar way
the B.C.s on the left end wall can be defined (Piechna [52]).

2.2. Determination of the Ports Timing
The procedure to determine the ports timing is as follows:
a) In the initialization step, the temperature and pressure in the channel are specified and the
velocity is set to 0.
b) At the beginning of the simulation, the outflow B.C. is applied to the outlet and the solid
wall B.C. is applied to the inlet.
c) Once the pressure at the inlet end of the channel is below the inlet pressure (p(1) < ptot,in),
the inlet B.C. is switched from the solid wall B.C. to the inflow B.C. and the Inlet
Opening time is recorded.
d) Once the condition of T(0.98N)new > 3T(0.98N)old is satisfied, the interface between the
mixture and the exhaust gas is considered as arriving at the outlet. Then the outlet B.C. is
switched from the outflow B.C. to the solid wall B.C. and the Outlet Closing time is
recorded.
25

e) Once the condition of u(0.02N)new < 0.2u(0.02N)old is satisfied, the hammer shock wave
is about to arrive at the inlet. Then the simulation stops and the Inlet Closing time is
recorded.
f) The Inlet Opening time, Outlet Closing time and Inlet Closing time multiplied by the
angular velocity yields the Inlet Opening angle AI_open, Outlet Closing angle OC_open
and Inlet Closing angle IC_closed.

2.3. Test Case and Results
A test case was tested to validate the code. The dimension of the computational domain and the
working condition is shown in Table 2-1, where ptot,in and Ttot,in are the total pressure and total
temperature at the inlet, p0 and T0 are the initial pressure and temperature in the channel, and
pout is the static pressure at the outlet.

26

Table 2-1: Dimension of the computational domain and the working condition of the test case for
the validation of the 1D code

IR (cm)

OR (cm)

ptot,in (Pa)

Ttot,in (K)

p0 (Pa)

T0 (K)

pout (Pa)

RPM

5

10

111325

300

709275

2100

101325

10000

The wave diagram was obtained from the TVDM scheme, shown in Figure 8. From the
pressure contour it can be clearly seen that after the opening of the outlet, an expansion wave
propagates towards the inlet, lowering down the pressure; and then it is reflected on the inlet wall
and the reflected expansion wave propagates towards the outlet; once arriving at the outlet, the
reflected expansion wave is reflected as a compression wave due to the subsonic outflow; during
the scavenging process expansion waves and compression waves travels back and forth; once the
outlet closes a hammer shock wave is generated and propagates towards the inlet. From the
temperature contour it can be seen that the expansion waves decrease the temperature and a clear
interface lies between the exhaust gas and the fresh mixture.
Table 2-2 lists the ports timing results from two schemes. The discrepancy shows that they
are consistent.

27

Table 2-2: Ports timing obtained by GT scheme and TVDM scheme
GT
-4

TVDM Discrepancy (%)

t

(×10 s)

1.3331

1.3340

0.07

t

(×10 s) 4.3046

4.3292

0.57

t

(×10 s) 6.0082

6.0538

0.76

IO

-4

OC
IC

-4

Analytical time for the head of the expansion wave to arrive at the inlet end is calculated as
follows:

t
=

L
=
a

L
=
γ RT

0.05
= 5.44 ×10−5 s
1.4 × 287 × 2100

Figure 9 shows the property profiles of two schemes at 5.44×10

-5

(2.31)

s. It can be seen that their

profiles overlap with each other, which shows a perfect consistency between them; the heads of
the expansion waves arrives at the inlet end, which is consistent with the analytic result.
Figure 10 shows the comparison between the results from GT scheme and TVDM scheme at
-5

7×10 s. The diagrams show that the results from two 1D codes match well.
Once the pressure at the inlet wall is just below the inlet pressure, the channel inlet opens
and the mixture is sucked into the channel. The mixture-exhaust interface can be obtained by
tracking the discontinuity location in the profiles. Figure 11 shows that the result from the GT
scheme and that from the TVDM scheme match well, but near the continuity the TVDM scheme
has some numerical oscillations.
Once the code sees the mixture-exhaust interface arrives at the outlet, the channel outlet

28

closes and a hammer shock is generated by bringing the mixture to rest. The location of the
shock wave can be obtained by tracking the discontinuity location in the profiles. Figure 12
shows that the wave front from the GT scheme is slightly ahead of that from the TVDM scheme,
which is because the outlet closed earlier in the GT scheme. When the inlet is closed, the
pressure in the channel is approximately 180000 Pa, which means the precompression ratio by
the hammer shock wave is 180000/111325 = 1.6.

2.1. Particle Tracking
It is useful to analyze the thermodynamic cycle of the engine through the T-s diagram and p-v
diagram. However, since the flow is unsteady, not all fluid particles experience the same
thermodynamic cycle and the flow field changes with time. Thus, the diagrams cannot be drawn
by acquiring data on various locations along the flow at one time instance as in a steady-state
flow. Tracking a particular particle is required to obtain the diagrams.
A numerical code was created to track the particles and draw the diagrams. In the tracking
process, the location of the fluid particle at the next time step is calculated based on the location
at the current time step and its current velocity, and the values of its properties at various time
step are recorded. At last, the diagrams are drawn by sequencing the property values in time.
Figure 13 shows the T-s and p-v diagrams obtained by tracking a single fluid particle from it
enters the channel to it leaves. The diagrams showed a clear Humphrey cycle, which is consistent
with the expectation.

29

CHAPTER 3. 2D SIMULATION
Since the 1D code does not take into account the 2D effect, such as the complex channel shape,
the leakage, gradual opening and closing of the ports and interaction between the channels, and
to apply more sophisticated combustion model and calculate the power and efficiency, the 2D
simulation is required for the design process.

3.1. Optimization Parameters
For the design in Figure 6, only the fresh air is drawn into the channels through the inlet. When
the channel closes, the fuel is directly injected to the channel and then ignited. However, the
mixing takes a lot of time and leaves short time for the combustion. Thus, to take advantage of
the fast speed of the premixed combustion, the air-fuel mixture is supplied through the inlet
which is called mixture inlet. In the configuration in Figure 6, the mixture will be ignited when
the channel with the combustion products opens to the mixture inlet and the hot gas contacts the
mixture. To prevent the backfire, another inlet is added before the mixture inlet. It is used for the
supply of the fresh air only, to create a buffer layer between the burned gas and the mixture. This
additional inlet is called air inlet. The configuration for the premixed combustion is shown in
Figure 14, along with the geometric parameters that need to be determined yet. The geometric
parameters contain the ports opening and closing angles, the inner and outer diameters of the
rotor, the inclination angles of the ports walls and the blade angles. The symbols and their
corresponding parameters are tabulated in Table 3-1.

30

Table 3-1: Symbols in Figure 14 and their meaning
AI_open

The opening angle of the air inlet

AI_closed

The closing angle of the air inlet

FI_open

The opening angle of the fuel inlet

FI_closed

The closing angle of the fuel inlet

EO_open The opening angle of the outlet
α_AI1

α_FI1

α_EO1

β_1

EO_closed The closing angle of the outlet

The right wall angle to the radial α_AI2

The left wall angle to the radial

direction of the air inlet

direction of the air inlet

The right wall angle to the radial α_FI2

The left wall angle to the radial

direction of the fuel inlet

direction of the fuel inlet

The right wall angle to the radial α_EO2

The left wall angle to the radial

direction of the outlet

direction of the outlet
β_2

Inlet blade angle

Outlet blade angle

The other parameters are the rotational speed, disc cycles, the channel number and the ratio
of the wall thickness to the channel width.
To create the channel wall in Gambit, a number of points on the wall are supposed to be
created first, which requires the knowledge of the coordinates of these points. For the first try, a
circular arc is chosen to be the channel shape (C1 shape). Its inlet blade angle is 100° and outlet
angle is 10°. Both walls have the identical shape. A Matlab code was made to calculate the
coordinates of the arc. Once the points are created, Gambit connects them using a NURBS curve.

3.2. Simulation methodology
For the simulations in this chapter, the ports timing is first obtained from the 1D code. Then the
2D simulations are performed using Fluent. Since the flow is unsteady, the transient solver is

31

used. Because the flow is mainly subsonic and the density-based solver is much slower than the
pressure-based one, the pressure-based solver is used. The standard k-ε model is used for the
turbulent flow. All walls are treated as adiabatic walls. Methane is chosen to be the fuel. The
Laminar Finite-Rate model is chosen for the combustion simulation, where the Arrhenius rate is
used for the reaction rate. Inlet Diffusion, Diffusion Energy Source, Full Multicomponent
Diffusion and Thermal Diffusion were turned on. The constant pressure specific heat is
calculated by the mixing law. The thermal conductivity is set to 0.0454 W/m-K and the viscosity
1.72×10

-5

kg/m-s. The mass diffusivity and the thermal diffusion coefficient follow the kinetic

theory. The spark ignition is simulated using a heater in the channel. In Figure 15 the rectangular
is the heater. Each channel has one heater near the inlet. In the beginning, the boundary type of
the rectangular is INTERIOR. When the channel rotates to a particular position, the heater is
turned on by switching its boundary type to RADIATOR and its temperature is set to 2000 K and
6

heat transfer coefficient is 10×10 W/m2-K. After 0.1 ms the heater is turned off by switching its
boundary type back to INTERIOR. The time step size is set to 5×10

-7

s. The working conditions

are listed in Table 3-2.

Table 3-2: Working conditions for the baseline design
Air inlet total pressure (barg)

0

Fuel inlet total temperature (K)

300

Fuel inlet total pressure (barg)

0.1

Outlet pressure (barg)

0

Air inlet total temperature (K)

300

RPM

10000

32

3.3. Optimization for a Given Rotor
After 4 rotations, the operation reaches the periodic mode. For a given rotor, the geometry of the
ports should be optimized for the best performance. The optimization is done in the following
aspects:
1) The air inlet opening angle should be adjusted to make sure when the channel inlet is
open, the pressure is below the ambient pressure, so that there is no backflow into the air inlet.
Figure 16(a) shows that when the channel is about to open to the air inlet, the pressure in the
channel is higher than that in the air inlet. As a results, when it is open to the air inlet, some
exhaust gas flows into the air inlet, as shown in Figure 16(b).
2) The wall angles of the ports should match with the blade angles to avoid the incidence
loss.
In Figure 17(a) the incoming flow impacts the suction surface of the blade, which generates
negative torque. The mismatching of the ports wall angles does not only decrease the torque
generation, but harms the fuel loading. In Figure 17(b) the mixture concentrates only on the
suction side of the blade so the channel cannot be fully filled. To prevent this, the inlet porting
walls should incline more.
The outlet walls should be diverging because when the channel is just open, the direction of
the absolute velocity at the channel outlet is opposite to the rotational direction and as channel
rotates the absolute velocity gradually turns to the rotational direction, as shown in Figure 18.
3) The air inlet is wide enough to generate an air layer between the hot gas and the fresh

33

mixture.
4) The closing angles of the outlet port and the fuel inlet are proper for the channel to be
filled up with the mixture without spill.
In Figure 17(b) the interface between the mixture and the air does not arrive at the outlet
end of the channel when the outlet is closed, so the outlet port and the fuel inlet port should be
enlarged.
5) No backfire should occur in the fuel inlet port.
Figure 19(a) shows that the hot burned gas is leaked into the following channels and the fuel
inlet. As a result, the mixture in the following channels and the fuel inlet is ignited by the leakage.
To prevent the backfire, the mixture needs to be ignited later so that the hot leakage can be
absorbed by the following channels.
6) The combustion should be complete before the channel is open to the outlet port.
Figure 20 shows that in this design when the channel is open to the outlet port, there is still
some unburned mixture left in the channel. To obtain a complete combustion, the spark location
should be moved to the center of the channel and one set of the ports should be removed to
ensure enough time for the combustion.

3.4. Baseline Design
After the geometry adjustment, the design was finalized. The final geometric parameters are
listed in Table 3-3.

34

Table 3-3: Final geometric parameters of the baseline WDE
ID (cm)

10

AI_open (°)

72

OD (cm)

20

AI_close (°)

112

Gap (cm)

0.04

FI_open (°)

113

α _AI1 (°)

50

FI_close (°)

170

α _AI2 (°)

50

EO_open (°)

25

α _FI1 (°)

50

EO_close (°)

150

α _FI2 (°)

50

β1 (°)

100

α _EO1 (°)

-50

β4 (°)

10

α _EO2 (°)

20

Channel number

30

Height (cm)

2.5

Channel-to-wall ratio

4

Figure 21(a) shows that the air inlet effectively created an air layer between the hot gas and
the fresh mixture. Figure 21(b) shows that the channel is filled up with the mixture when the
channel is closed and the combustion is complete before the channel is open. Figure 21(c) shows
that the relative velocity of the incoming flow follows the blade shape.
Table 3-4 shows the power and efficiency of the baseline engine in the second rotation. The
power is calculated from the equation below
P = τ meanω

(3.1)

where τ mean is the mean torque during the chosen rotation and ω the angular velocity of the
rotor. The efficiency is calculated from the equation below

η=

P
m mixtureYCH4 LHVCH4

35

(3.2)

where m mixture is the average mass flow rate through the fuel inlet, YCH4 is the mass fraction
of methane in the mixture and LHVCH4 is the lower heating value of methane. Theoretically, the
efficiency of Humphrey cycle without a pre-compression ratio of 1 is about 30% which is much
higher than the simulation. Thus, the energy loss should be examined to enhance the efficiency.

Table 3-4: Power and efficiency of the baseline design
Cycle Power (kW) Efficiency (%)
2

2.3712

1.23

Figure 22 shows the power history in the second rotation. Due to the unsteady nature, the
power fluctuates with time.

36

CHAPTER 4. METHODS TO IMPROVE THE EFFICIENCY
4.1. Investigation into the Energy Usage
4.1.1.

1D Investigation

First, the particle tracking technique based on the 1D simulation was used to investigate the
thermodynamic process that different fluid particles experience.
Figure 23 shows the T-s diagrams of 10 fluid particles that enter the channel at different
time. From the figure it can be found that due to the unsteady nature of the flow, not all the
exhaust gas particles leaves the channel at the ambient pressure, meaning that not all gas
experiences a full Humphrey cycle. The gas that leaves first is at a high pressure and temperature,
which is one of the sources of the energy lost and a low efficiency.

4.1.2.

2D Investigation

To verify the conclusion from the 1D investigation and take into account the 2D effect (channel
curvature, gradual opening and closing of the channel), a further investigation was performed by
2D simulation.
Figure 24 shows the geometric configuration of the 2D simulation. To shorten the
simulation time, only one channel is examined. The rotor has two sets of ports, so the channel
experiences two cycles in each revolution. To eliminate the influence of the leakage, the gaps
between the rotor and the stator were removed.
In this investigation only the thermodynamic characteristic is of interest and, a lot of

37

information from the combustion model is redundant, such as the flame propagation and the
concentration of the species. Thus, a simpler model for the combustion can be adopted to
accelerate the simulation. The species transport and reaction model was turned off and the air
was chosen as the working material. Since the gradualness of the heat release by the combustion
does not affect the power generation in this perfectly closed channel, the heat release can be
modeled as a sudden rise in temperature and pressure. This was done by patching a high
temperature and high pressure to the channel when it is fully closed. In this investigation, 2100 K
was chosen as the patching temperature assuming that after the combustion the temperature in
the channel increases to 2100 K. if the pressure is patched the arbitrarily, the density in the
channel would change and consequently the mass is not conserved. To satisfy the mass
conservation, the channel pressure should be patched according to the ideal gas law:
p3 =

2100K
p2
T2

(4.1)

where T2 and p2 are the temperature and pressure before the patching and p3 is the patching
pressure.
As the thermodynamic analysis in section 1.3, the channel is taken as the control volume
and the integral form of the energy equation Eq. (1.1) is applied.
The simulation ran for 6 working cycles (3 revolutions). The results are shown in Table 4-1.
The net work done by the fluid to the channel was calculated by the following equation

Wnet = τ meanωTcycle
The heat addition was obtained by

38

(4.2)

Q
= E3 − E2

(4.3)

where E2 is the internal energy in the channel before the patching and E3 is that after the
patching.
The cycle efficiency is

η=

Wnet
Q

(4.4)

The kinetic energy is
=
k H tot − H

(4.5)

The potential efficiency is

η potential =

Wnet + ( kout − kin )
Q

(4.6)

The discrepancy of the energy equation is

=
∆ energy

(

)

Q − H tot ,out − H tot ,in + Wnet
×100%
Q

(4.7)

The discrepancy of the mass conservation equation is

=
∆ mass

m in, mean − m out , mean
m in, mean

39

×100%

(4.8)

Table 4-1: Results of the 2D simulation for the analysis of the energy loss
Cycle

1

2

3

4

5

6

Wnet (J)

0.8423

0.9905

0.9885

1.0014

1.0061

1.0061

𝑸 (J)

23.9362

26.6363

26.3928

26.7855

26.8005

26.8326

𝜼 (%)

3.52

3.72

3.75

3.74

3.75

3.75

Htot,in (J)

0.2161

0.2115

0.1036

0.1037

0.1002

0.1000

Hin (J)

0.1089

0.0965

-0.0268

-0.0271

-0.0307

-0.0305

kin (J)

0.1072

0.115

0.1304

0.1308

0.1309

0.1305

Htot,out (J)

-22.7106

-25.463

-25.2206

-25.5301

-25.5395

-25.5518

Hout (J)

-18.4599

-20.4108

-20.2095

-20.4261

-20.4228

-20.4309

Kout (J)

-4.2507

-5.0522

-5.0111

-5.1040

-5.1167

-5.1209

ηpotential (%)

20.83

22.25

22.24

22.31

22.36

22.35

Δenergy (%)

2.50

1.48

1.09

1.34

1.32

1.40

m in,mean (g/s)

8.269

8.3601

8.4485

8.4705

8.4127

8.3443

m out,mean (g/s)

-7.4641

-8.4208

-8.3621

-8.4647

-8.4032

-8.3643

Δmass (%)

9.73

-0.726

1.02

0.0685

0.113

-0.240

pcompression (Pa)

125400.7

124597.9

124567

124861

125140.6

125068.3

pcompression / pin

1.126

1.119

1.119

1.122

1.124

1.123

For the sixth cycles, the discrepancy of the energy equation is 1.40% and that of the mass
40

conservation equation is -0.24%, which means that both of the energy and mass accumulated in
the channel can be negligible and the working mode of the sixth cycle can be considered as a
periodic mode.
In section 1.3, to achieve the efficiency of the Humphrey cycle, it was assumed that the
incoming kinetic energy is equal to the outgoing one. However, in this simulation, the latter
(5.1209 J) is much greater than the former (0.1305 J), meaning that a lot of energy was released
to the ambience in the form of kinetic energy without being utilized. If the outgoing kinetic
energy could be reused, the efficiency would increase from 3.75 % to 22.35 %.
Moreover, in the sixth cycle after the patching, the mean pressure in the channel is
805654.83 Pa. In the ideal Humphrey cycle, when the gas expanded isentropically from p3 =
805654.83 Pa and T3 = 2100 K to p4 = 1 atm, the temperature should be lowered down to T4 =
1162.64 K calculated from the isentropic relationship
γ −1
 p4  γ

T4
=

T3  p3 

(4.9)

Thus, the corresponding specific enthalpy is

(

)

h4 = c p T4 − Tref = 1006.43 × (1162.64 − 298.15 ) = 869699 ( J/kg )

(4.10)

From Figure 25 it can be seen that the specific enthalpy at the outlet of the channel
gradually drops down after the channel is open, and it is above the theoretical value in the first
half of the expansion. Theoretically, during the expansion process (0.0129435 s to 0.0133225 s),
the outgoing enthalpy should have been

41

 4 ∆=
H=
t 0.0365999 × ( 0.0133225 − 0.0129435 ) × 869699
= 13.8071( J )
4 mh

(4.11)

However, from the simulation, this value is 16.3922 J which is higher than the theoretical
one, meaning that not all the gas is fully expanded. In fact, from the Figure 25 it can be seen that
the gas leaving first is under-expanded and that leaving later is over-expanded. The enthalpy
analysis shows that the gas is under-expanded in terms of the net effect.
Based on Eq. (4.9) for the Humphrey cycle with the precompression ratio of 1 (Lenoir
cycle), if the temperature and pressure after the heat addition are T3 = 2100 K and p3 = 7 atm, the
temperature after the expansion is T4 = 1205.7 K. The thermal efficiency of an air-standard can
be calculated as follows:
1 γ
η =−

T4 − T1
1205.7 − 300
1 1.4 ×
=−
=29.6%
T3 − T2
2100 − 300

(4.12)

In sum, there are two sources of the energy loss in this simulation. One is the difference
between the outgoing kinetic energy and the incoming one, and the other is the under-expansion
of the gas. Thus, the method to improve the efficiency would be further expanding the gas in the
channel, reducing the outgoing kinetic energy, increasing the incoming kinetic energy and
reusing the energy from the exhaust gas.

4.2. Investigation into the Channel Shape
To improve the efficiency, more energy should be extracted from the rotor. To exploit the
potential of energy extraction from the channel, more advanced channel shapes need to be
explored.

42

4.2.1.

Torque Generation Mechanism for Unsteady Flow

Although researchers intended to use the curved channel in wave engines to generate power, the
principle of the power generation has not been explained clearly. Since the WDE is an unsteady
flow device, the power generation principle should be distinguished from the traditional steady
flow turbine.
For a steady flow turbomachinery, when a fluid particle passes through a blade passage, it is
subject to the centripetal force provided by the pressure gradient pointing to the center of blade
curvature. The pressure gradient increases with the curvature, the velocity and the density. One
way to increase the pressure gradient is to increase the curvature. Based on this idea the C1
channel is modified to the C2 shape to obtain a higher efficiency. The curve of C2 is obtained by
increasing the inlet blade angle from 100° to 150° while keeping the outlet blade angle of 10°.
Both walls are of the identical shape and the volume of the C2 channel is identical to that of the
C1 channel. The geometry of the C2 channel is shown in Figure 26.
To compare their performance in a steady flow, two steady state simulations were done for
them respectively. In the simulations, the gauge total pressure at the inlet is set to 0.1 bar and that
at the outlet is 0.
Figure 27 shows the pressure contour in the tested channels. Due to the curvature, the
pressure on the pressure side of the channel is higher than that on the suction side in both cases.
For the C1 channel the flow generates 7.09 N·m of torque, and the torque on the C2 channel is
14.59 N·m. The torque of C2 channel is approximately twice as much as the C1 channel, which
is consistent with the principle above. However, since the WDE operates in an unsteady flow, the
43

principle for the steady flow blade design might not be valid. Thus, the performance of power
generation in unsteady flow must be examined.
For a steady flow turbine, the torque on the blades does not vary with time. In contrast, the
torque on the channel walls of a WDE fluctuates a lot. The moment coefficient history in the
sixth cycle of the energy loss simulation is shown in Figure 28, where the moment coefficient

Cm =

τ

(4.12)

1
2
ρ ref vref
Aref Lref
2

is a non-dimensional torque.
In the torque history diagram, the area between the curve and the time axis shows the
energy extracted from the flow. It can be seen that the major extraction is done in the expansion
process. The energy extracted in each process is tabulated in Table 4-2. From the table it is
noticed that the major power is generated in the expansion process when the outlet is open and
the inlet is closed. The power generated in the other processes is negligible. Thus, in this section
the energy extraction performance of various channel shapes in unsteady flow is investigated
only for the expansion process.

Table 4-2: Work generated in each process in the sixth cycle
Expansion

Scavenging

Compression

Combustion

Starting time (s)

0.0159435

0.016322

0.0177695

-

Ending time (s)

0.016322

0.0177695

0.017915

-

Work generated (J)

1.051

-0.02726

-0.02078

0.002889

Contribution to total work (%)

104.5

-2.709

-2.065

0.2872

44

Several channel shapes were tested. In each case, the computation domain consists of a
rotating channel and a stationary outlet port. In the beginning of the simulation the channel was
fully closed with the air at 1 atm and 300K. Then the channel temperature is patched to 2100 K
and the pressure is patched based on Eq. (4.1). The simulation proceeds as the channel rotates at
10,000 RPM. Finally, the simulation stops when the pressure at the inlet wall of the channel
decreases to the ambient pressure.
To compare various channel shapes, the C1 channel was chosen as the baseline shape. The
geometry of the C1 channel is shown in Figure 29.
The work done by the flow to the channel wall is calculated using Eq. (4.2). The heat
addition is the internal energy increase due to the patching. The efficiency is calculated by
dividing the energy generated by the heat addition. The tests for the C1 channel and C2 channel
were done, and the results are shown in Table 4-3. However, the efficiency of C2 is not higher
than the C1 channel as in the steady flow tests. Thus, the unsteady factors must play an important
role in the energy extraction and the traditional blade design for steady flow does not apply.

Table 4-3: Results of the unsteady simulations for C1 and C2 channel
Shape

C1

C2

Duration (ms) 0.3835 0.3445
3.38

Efficiency (%)

3.13

To obtain a guidance to improve the channel shape, the power generation principle should

45

be investigated first.
The source of the torque on the channel wall is the pressure and the viscous force exerted by
the fluid. However, compared to the pressure, the effect of the viscous force is negligible. Thus,
it is only necessary to investigate the pressure in the channel. The pressure on the pressure side
generates positive torque, while the pressure on the suction side generates negative torque. So the
net torque depends on the pressure differential.
Previous research has shown that the front of the acoustic/shock wave propagating through
a curved channel is approximately perpendicular to the channel wall.
The similar phenomenon has been observed in the C1 channel. As shown in Figure 31, as
the channel gradually opens, the expansion wave originates at the trailing edge of the upper wall,
and then the wave front adjusts to be perpendicular to the channel walls.
During propagating towards the inlet, the wave front is always perpendicular to the walls,
shown in Figure 32.
While the reflecting wave interferes with the incident wave, the front is convex to the inlet
wall, shown in Figure 33.
The torque on the channel walls is generated by the pressure distribution. For simplicity, a
rotating straight channel was tested to explain the principle of torque generation by the expansion
wave. In the beginning of the simulation, the channel contains the gas at 7 bar and 2100 K, and
then the outlet suddenly opens to the ambience. Figure 34 shows the pressure contours during the
expansion wave travels towards the inlet end of the channel.
Figure 34 shows the expansion wave propagates towards the inlet end of the channel,
46

causing a pressure gradient after its head. In this analysis, a positive pressure gradient is defined
as the pressure decrease in the direction from the inlet to the outlet. If the pressure decreases in
the opposite direction, the pressure gradient is negative. It can be found that the isobaric lines are
perpendicular to the channel walls. Two points on the channel walls, point a and point b, are used
for analysis. They are the intersections of the circumferential line and the two walls, i.e., they are
at the same radius. Since the walls are not radial, each isobaric line has an angle with respect to
the circumferential direction. Due to this angle and the pressure distribution, there is a pressure
difference between point a and point b. In this case where the inlet angle of the channel wall is
smaller than 90°, the pressure at point a is higher than that at point b. In the same way, all points
on the upper wall are at higher pressure than their corresponding points on the bottom wall,
which generates positive net torque. If the blade angle is greater than 90°, the net torque would
be negative. A more general principle of the torque generation can be expressed by Table 4-4.

Table 4-4: General principle of torque generation by travelling expansion wave
Blade angle < 90° Blade angle > 90°
+ pressure gradient

+

-

- pressure gradient

-

+

The result that the C2 channel has a lower efficiency than C1 channel can be explained by
this theory. Figure 35 shows the pressure contours of the C2 channels. It can be observed that
near the outlet the blade angle is smaller than 90° and near the inlet the blade angle is greater

47

than 90°. Consequently, the torque generated by the part near the outlet is positive, while that
near the inlet is negative. According to the positive efficiency, the net torque it positive. Near the
inlet, the blade angle of the C2 channel is greater than that of the C1 channel, so the angle
between the isobaric line near the inlet and the circumferential direction is greater in the C2
channel, which causes greater negative torque than the C1 channel.
To further justify this principle, two new channel shapes were proposed. The first one is a C
shape channel with 70° of the inlet blade angle and 10° of the outlet one, which is shown in
Figure 36.
The second one consists of two parts - a straight line and a circular arc, and the straight line
is tangent to the arc. This shape is named J shape. Figure 37 shows two channels of the J shape.
Both of them have the inlet blade angle of 100° and outlet blade angle of 10°. For J1 channel, the
radius of the tangent point of the straight portion and the arc is 7 cm. for J2 channel, the tangent
point is at 8 cm of radius. If the length of the straight portion reduces to 0 cm, the J shape
becomes the C shape.
The results of these three channels are tabulated in Table 4-5. Among C1, J1 and J2
channels, J2 has the lowest efficiency while C1 has the highest one because J2 has the longest
straight portion where the negative torque is generated, while the C1 has the smallest. Since the
blade angle is smaller than 90° everywhere for C3 channel, the torque is always positive. Thus,
C3 has the highest efficiency among all channels that have been tested so far.

48

Table 4-5: Results of the single channel simulations for C3, J1 and J2 channels
Shape

C3

J1

J2

Duration (ms) 0.4515 0.3290 0.2935
Efficiency (%)

4.2.2.

3.89

2.94

2.33

Effect of the Outlet Opening

If we plot the relationship between the efficiency and the expansion duration for the tests in the
previous section (see Figure 38 ), it can be found that the efficiency increases with the expansion
duration.
One important geometric parameter that influences the expansion duration is the outlet
opening, which is defined as the smallest distance from the trailing edge of the upper wall to the
surface of the lower wall, shown in Figure 39.
From the relationship between the expansion duration and the outlet opening (Figure 40), it
can be found that the expansion process is slowed down by the smaller outlet opening.
Combining with Figure 38, the smaller outlet opening gives higher efficiency. Based on this idea,
the J2 and C3 channels are modified to have smaller outlet openings.
The modification method is to rearrange the curvature of the upper wall while maintaining
the outlet blade angle. The geometries of the modified channels (MJ2 and MC3) are shown in
Figure 41. Their outlet openings, expansion durations and efficiencies are listed in Table 4-6.
From the results, it can be seen that both of them have higher efficiency than the channel tested
in the previous section and their outlet openings are smaller than those channels, which justifies
49

the relationship between the outlet opening and the efficiency. However, even if the MC3
channel has a bigger outlet opening and shorter expansion duration than the MJ2 channel, its
efficiency is still higher, which can be explained by the torque generation principle describe in
the previous section.

Table 4-6: Outlet openings, expansion durations and efficiencies of the single channel
simulations for J2, C3, MJ2 and MC3 channels
Shape

J2

MJ2

C3

MC3

Outlet opening (cm)

0.7971 0.1500 0.4350 0.1515

Expansion duration (ms) 0.2935 1.5947 0.4515 1.275
Efficiency (%)

2.33

4.77

3.89

5.12

For all channels that have been tested so far, both walls have the identical inlet blade angle
and the identical outlet blade angle. A convergent channel with two different shape walls was
tested additionally. The geometry is shown in Figure 42, and the blade angles are listed in Table
4-7.

Table 4-7: Blade angles of the convergent channel
Upper wall Lower wall
Inlet angle (°)

75

90

Outlet angle (°)

10

10

Its outlet opening, expansion duration and efficiency are listed in Table 4-8. Due to the
50

bigger outlet opening, the expansion duration of the convergent channel is shorter than those of
MJ2 and MC3. However, the efficiency is higher than MJ2, which is because the inlet blade
angles are smaller than those of MJ2. The efficiency is lower than MC3 is resulted from that the
blade angles and the outlet opening are bigger than those of MC3.

Table 4-8: Results of the single channel simulations for the convergent channel
Outlet opening (cm) Expansion duration (ms) Efficiency (%)
0.2338

0.7860

4.82

Comparing the torque histories of all tested channels (see Figure 43), it can be found that
for the short expansion the torque has only one peak, while for the long expansion the torque has
multiple peaks. As the expansion gets longer, the torque oscillation becomes clearer.
The torque oscillation can be explained by the variation of the pressure gradient with time.
Take the MC3 channel for example. When the channel outlet opens, an incident expansion wave
propagates into the channel, causing the positive pressure gradient behind the wave front and
accelerating the mass flow out of the channel. As the incident expansion wave propagates
towards the inlet end, the area with the positive pressure gradient grows (Figure ). Since the
blade angle is smaller than 90°, according to the torque generation mechanism, the positive
torque increases in this duration. The torque can be illustrated directly by the profiles of pressure
as a function of radius on the pressure surface and the suction surface (Figure ). In the profile

51

plots, the area under each curve is the tangential component of the force on the corresponding
surface and the area between the black curve and the red one gives the net tangential force on the
channel (shaded area). The force is positive in the regions where the black curve is above the red
curve and negative where the black curve is below the red one. Figure shows that the region of
positive force grows as the incident expansion wave propagates. When the head of the expansion
wave arrives at the inlet end, the positive torque is generated in the entire channel and reaches
the first peak. Therefore, the propagation of the incident expansion wave is the main reason for
the first rise of the torque.
Besides the propagation of the incident expansion wave, the gradual opening is another
reason. Figure 45 shows that there is a relatively high pressure region (the circled area in the
pressure contour) near the trailing edge of the suction surface, and the gauge pressure decreases
to the ambient pressure at r = 0.099 cm and then increases to 2.5 bar, which is because the
bottom corner of the channel outlet is still closed. Compared with the imaginary scenario where
the pressure remains ambient pressure beyond r = 0.099 cm, the net force is smaller (the shaded
area in the bottom diagram) in the real case and hence a smaller positive torque. This is the
undesirable effect of the gradual opening. To reduce this effect, the opening time should be
shortened by either increasing the blade tip speed or decreasing the channel circumferential
opening (the circumferential length of the channel outlet).
Figure

shows that when the channel outlet is fully open, the high pressure region

disappears and the gauge pressure near the trailing edge of the suction surface decreases to -1 bar.
From Figure one can see that the flow is supersonic outside the channel opening. When the
52

supersonic flow sees the trailing edge of the channel bent away from the flow, a Prandtl-Meyer
expansion wave is generated on the bent edge further decreasing the pressure. Compared with
the scenario where the pressure remained atmospheric in the region where the Prandtl-Meyer
expansion wave is located, the Prandtl-Meyer expansion wave generates additional force (shaded
area in the middle diagram) and hence an additional positive torque. This is the benefit of the
supersonic outflow and the curved trailing edge. Thus, the shrinking of the high pressure region
as the channel outlet opens is the minor reason for the first rise of torque.
Once reaching the inlet end, the incident expansion wave reflects. The reflected expansion
wave further decreases the pressure behind its head and decelerates the mass flow. As
propagating towards the outlet, the reflected wave interferes with the rest of the incident
expansion wave, forming a nonsimple region. In this case and other the long expansion cases,
such as the convergent channel and MJ2, the pressure gradient is negative (Figure 47) in the
nonsimple region, and the peak pressure occurs at the interface between the simple and
nonsimple regions. According to the torque generation principle, the negative pressure gradient
with the blade angle smaller than 90° generates negative torque. In the short expansion cases,
although the negative pressure gradient is not observed, the positive pressure gradient in the
nonsimple region is very small. As the reflected expansion wave propagates towards the outlet,
the nonsimple region grows. Thus, during this process the net torque on the channel walls
decreases. Due to the fact that the flow is choked at the outlet, the tail of the incident expansion
wave cannot propagate into the channel. Therefore, when the head of the reflected wave arrives
at the outlet, the nonsimple region is all over the channel. However, near the outlet of the channel,
53

the positive pressure gradient caused by the incident expansion wave is greater than the negative
pressure gradient caused by the reflected expansion wave, the net pressure gradient is still
positive. Therefore, the net torque on the walls is still positive.
As the reflected expansion wave propagates, the mass velocity in the nonsimple region
decreases to zero and then to negative values due to the inertia. In other word, the gas flows
backward (Figure 48). When the backflow impacts on the inlet end wall, a compression wave is
generated at the inlet end. The compression wave and the incident expansion wave cause a
positive pressure gradient behind the compression wave, shown in Figure 49. As the compression
wave propagates towards the outlet, the positive pressure gradient zone grows, and consequently
the positive torque increases near the inlet end. Thus, the net torque in the whole channel
increases. Eventually, the pressure gradient becomes positive everywhere in the channel and
torque oscillation cycle goes back to the state when the incident expansion wave propagates
towards the inlet.
The process described above repeats several times, and the positive and negative pressure
gradients occur near the inlet of the channel alternately, causing the torque oscillation. The
torque oscillation phenomenon is very obvious in MJ2, MC3 and the convergent channels; in the
case of C3 channel, the torque oscillates less than the three long expansion cases but more than
the other short expansion cases. Although the smaller channel opening prolongs the expansion
process, the reflected expansion wave has a negative effect on the torque generation.
Another conclusion can be drawn by the observation of the pressure profiles. In all
processes, the maximum tangential force (the area between the black curve and the red one in
54

Figure , Figure 47 and Figure 49) occurs near the outlet of the channel. This force generates
greater torque when the channel outlet is at the outer radius than at the inner radius due to a
longer force arm at the outer radius. Thus, the channel outlet should be at the outer radius so that
the force can be made used of more efficiently.

4.3. Implementation of the New Shape
Although the longer expansion gives the higher efficiency, it has some disadvantages. For
example, at a certain rotational speed, the total time for the expansion, scavenging, compression
and combustion is fixed. The long expansion channels have smaller outlet openings, which
causes a longer scavenging. Since both the expansion and scavenging times are longer, the time
left for the combustion is shorter and the combustion may not be complete. In this section, the
convergent channel is chosen to be implemented on a WDE, for it has a high efficiency and a
relatively shorter expansion.
To compare with the convergent channel engine, a baseline simulation was performed first.
The geometry of the baseline engine was taken to be tested. However, in this simulation the
combustion behavior is not of interest, but the aerodynamic performance of the engine. Thus, the
species transport and reaction model were turned off and the patching method was used for
combustion modeling as in section 4.1.2. To reduce the influence of leakage, the channel is
patched just before it opens to the outlet port as shown in Figure 50.
The histories of the torque of the first three cycles (Figure ) show that the third cycle reaches
the periodic operation mode. Table 4-9 lists the power and efficiency in the third cycle.
55

Table 4-9: Results of the simulation with patching method for the baseline design
Cycle Power (kW) Efficiency (%)
3

1.6476

2.21

The geometry of the WDE with the convergent channels is shown in Figure 52. The
histories of torque (Figure 53) show that the second cycle reaches the periodic mode.
Comparison between Table 4-9 and Table 4-10 shows that the efficiency increases by 81%.

Table 4-10: Results of the WDE with the convergent channels
Cycle Power (kW) Efficiency (%)
2

2.4803

3.96

4.4. Effect of the Peak Pressure
4.4.1.

1D Investigation

To investigate the reason for the under-expansion of the exhaust, two cases were tested using the
1D code. The first one is the test case in section 2.3.
Figure 54 shows that the outflow is choked, so the tail of the expansion wave cannot
propagate into the channel. The flow is under-expanded.
In the second case, the initial pressure is reduced to 3 bar. Figure 55 shows that in this case

56

the outflow became subsonic, and the entire expansion wave was confined in the channel. Thus,
reducing the pressure in the channel after the combustion leads to a way to fully expand the gas.

4.4.2.

2D Investigation

To justify the result from the 1D simulation above, two 2D cases were tested using the geometry
of C shape 3. Their peak pressures were 5 atm and 4 atm respectively. Table 4-11 lists the results.
The results show that the efficiency increases as the peak pressure decreases, which is consistent
with the 1D simulation. However, as the peak pressure decreases, the power decreases as well.
Thus, the determination of the peak pressure depends on the priority of the power and efficiency.

Table 4-11: Results of the initial pressure tests

4.4.3.

Peak pressure (atm)

7

5

4

Power (kW)

2.0626

1.7403

1.4567

Efficiency (%)

3.89

4.66

5.08

Full Engine Test

C2 channel was used on the full engine test. Figure 56 shows the geometry of the engine. In this
simulation the patching temperature is 2100 K and the pressure was patched accordingly. The
results of the second rotation are shown in Table 4-12.

57

Table 4-12: Results of the simulation with 2100 K patching temperature for the WDE with C3
channels
Cycle Power (kW) Efficiency (%)
2

2.9905

1.94

Then the simulation was re-run with 1200 K of patching temperature and the results the
second rotation are shown in Table 4-13. Compared to the previous simulation, the case with a
lower pressure and temperature after combustion has a higher efficiency, which is consistent with
the 1D and 2D investigation. In practice, the lower combustion pressure can be achieved by
filling a lean mixture or partially filling the channel.

Table 4-13: Results of the simulation with 1200 K patching temperature for the WDE with C3
channels
Cycle Power (kW) Efficiency (%)
2

1.7735

2.56

4.5. The Effect of the Blade Tip Speed
4.5.1.

Single Channel Investigation

According to the energy loss analysis, the most important loss is the significant outgoing kinetic
energy. If the rotor is stationary, when the channel opens to the ambience, the burned gas
exhausts at a very high velocity due to the pressure difference and all energy is released to the

58

ambience without being used. As the rotor accelerates, the blade tip speed increases. The relative
velocity of the exhaust gas mainly depends on the pressure ratio of the gas and the ambience. In

the case shown in Figure 57, under the condition that the tangential relative velocity wθ is to

�⃑ and the magnitude of the tip speed U is smaller than that of the
the opposite of the tip speed 𝑈
tangential relative velocity wθ, the magnitude of the tangential absolute velocity cθ decreases as

the tip speed increases. Since the radial absolute velocity cr is the radial relative velocity wr
which does not change with the tip speed, the magnitude of the absolute velocity c decreases as
the tip speed increases. Therefore, the outgoing kinetic energy decreases. For the expansion
process when the inlet is closed, there is no energy coming in, and the outgoing enthalpy depends
only on the pressure ratio and does not change with the tip speed, so a smaller outgoing kinetic
energy means a higher energy extraction and hence a higher efficiency. However, if the
magnitude of the tip speed is greater than that of the tangential relative velocity, the tangential
absolute velocity increases with the tip speed and hence an increasing kinetic energy. Thus, the
performance of the energy extraction increases up to a maximum and then decreases as the tip
speed increases.
Besides, as discussed in section 0, increasing the blade tip speed reduces the undesirable
effect of the gradual opening and benefits the torque generation.
There are two ways to increase the tip speed – increasing the outer diameter and increasing
the rotational speed.
A series of single channel simulations were performed to justify the relationship between
the efficiency and the rotor dimension. The convergent channel was taken for the test. In this test,
59

the radial distance between the outlet and the inlet remained 5 cm in each case, and three outer
radii were tested – 10 cm, 15 cm and 20 cm. Table 4-14 and Figure 58 show that both power and
efficiency increase with the channel diameter within the tested range.

Table 4-14: Results of the rotor dimension tests
OR (cm)

10

15

20

IR (cm)

5

10

15

Power (kW)
Efficiency (%)

1.3506 2.1152 2.7754
4.82

7.14

10.32

The C3 channel was taken to examine the relationship between the efficiency and the
rotational speed. In this test, five speeds - 1000 RPM, 5000 RPM, 10000RPM, 15000 RPM and
20000 RPM - were tested. The initial temperature and pressure are identical to those in section
4.2. Their results of five rotational speeds are listed in Table 4-15. Figure 59 shows that both
power and efficiency increase with the rotational speed within the tested range. Another
advantage of high rotational speed is that the opening and closing of the channel is faster, which
is favorable to the generation of expansion wave and hammer shock wave.

60

Table 4-15: Results of the rotational speed tests
1000

5000

10000

15000

20000

Duration (ms)

0.7305

0.479

0.4515

0.4416

0.4365

Power (kW)

0.0307

0.7642

2.0626

3.5081

4.7749

Efficiency (%)

0.0937

1.53

3.89

6.47

8.71

RPM

To adopt a higher rotational speed the channel dimension should be adapted accordingly.
The higher rotational speed shortens the time for combustion, which is unfavorable for the
completeness of the combustion. For example, if the rotational speed is doubled, the optimistic
case is that all times for expansion, scavenging, compression and combustion are halved.
However, this is not the real case. Table 4-15 shows that the expansion duration decreases as the
rotational speed increases due to the increasing centrifugal force, but the variation is insignificant.
From 5000 RPM to 10000 RPM, the rotational speed increases 4 times but the duration decreases
only 8.87%. The duration of 1000 RPM is irregularly long because of the gradualness of the
outlet opening. Thus, it is reasonable to assume that the wave time (the expansion and
compression durations) does not change with the rotational speed. To examine the variation of
the scavenging time with the rotational speed, a 1D simulation was performed. The geometry and
working conditions are the same as the test case shown in Table 2-1, only with the rotational
speed changed to 20000 RPM. From Table 4-16, the scavenging time can be calculated. The time
-4

for 10000 RPM is 4.3292 – 1.3340 = 2.9952 (×10 s), and that for 20000 RPM is 4.0657 –
-4

-4

1.3111 = 2.7546 (×10 s) which is 8.03 % less than 2.9952 (×10 s). Consequently, for a doubled

61

rotational speed the time for combustion is less than half of the original time. To solve this
problem, the channel length should be halved so that the wave time and scavenging time can be
halved and the time required by complete combustion is halved as well. Therefore, for a given
blade tip speed, a rotor with the original outer diameter but half length and double rotational
speed can use the same ports timing as the rotor with the original length and rotational speed but
a doubled outer diameter. In other words, the dimension of the engine can be scaled without
changing the ports timing and the rotational speed should be scaled accordingly. The effect of the
scaling on the efficiency will be examined in section 4.6.4.

Table 4-16: Comparison of the ports timing of different rotational speeds

RPM

4.5.2.

-4

-4

-4

tIO(×10 s) tOC(×10 s) tIC (×10 s)

10000

1.3340

4.3292

6.0538

20000

1.3111

4.0657

5.8655

Full Engine Investigation

To examine if this principle applies for the full engine, a WDE with 20 cm ID and 30 cm OD was
created. Figure 60 shows its geometry. Since its perimeter is longer, it can accommodate more
channels. In this design there are 30 channels on the rotor.
Comparison between Table 4-17 and Table 4-10 shows clearly that the efficiency increases
by 26.5% by enlarging the rotor.

62

Table 4-17: Results of the WDE with 30 cm OD
Cycle Power (kW) Efficiency (%)
2

2.8632

5.01

4.6. Reuse of the Energy from the Exhaust Gas
4.6.1.

Return Channel

Since a significant amount of energy in the exhaust gas is not fully used, the efficiency of the
baseline design is much lower than the theoretical value of Humphrey cycle. If the exhaust gas
can be inducted back to the channel, some energy from the gas can generate additional power.
Adding a return channel serves this purpose.
Figure 61 shows the configuration of a WDE with a return channel. Table 4-18 lists its
power and efficiency. Comparing Table 4-10 with Table 4-18, the results show that the efficiency
increases by 56.1% by adding a return channel.

Table 4-18: Results of the WDE with a return channel
Cycle Power (kW) Efficiency (%)
2

3.9943

6.18

63

4.6.2.

Two-Stage Rotor

When the gas passes through the long return channel, pressure loss occurs in this process. It
would be beneficial if the exhaust gas can generate power immediately after leaving the channel.
Like a conventional gas turbine, adding a second rotor outside the engine can meet this
requirement.
Figure 62 shows the configuration of a two-stage engine. In this design several guide vanes
were added to the exit port of the first stage (inner rotor) to deflect the exhaust gas to the second
stage (outer rotor). The two rotors are connected to the same shaft. The combustion takes place
in the first stage and the second stage is used for energy extraction only.
When the inlet of the turbine blade passage opens to the outlet of the guide vane passage, a
compression wave propagates into the turbine passage. Thus, the torque generation mechanism
of the external turbine is similar to that of the inner rotor, and 70° was used for the inlet blade
angle and 10° for the outlet angle. The power and efficiency of the second rotation are listed in
Table 4-19. Comparison between Table 4-10 and Table 4-19 shows that the efficiency increases
by 106.8%. Comparing Table 4-18 with Table 4-19, one can find that the efficiency increases by
32.5%. Thus, using an external turbine is a more effective way to improve the efficiency.

Table 4-19: Results of the two-stage engine
Cycle Power (kW) Efficiency (%)
2

5.2240

8.19

64

4.6.3.

Configuration of the Guide Vanes and the Turbine Blades

Since the two-stage engine has the highest efficiency so far, more attention is paid to the
development of this configuration.
Figure 63 shows that when the blade passage (passage A) first opens to the guide vane
passage the compression wave propagates into the blade passage and the pressure gradient is
positive. However, as the blade passage (passage B) rotates, its inlet is exposed to a lower
pressure region and an expansion wave propagates into the blade passage causing a negative
pressure gradient near the inlet generating negative torque. As the blade rotates, the pressure
gradient near the inlet turns positive and negative alternately. Thus, the inlet angle smaller than
90° is not always beneficial. Since the blade passages have through-flow in them when they
generate torque, the operation of the blades is closer to that of a conventional turbine. Therefore,
the conventional turbine geometry may help improve the efficiency.
Figure 64 shows that due to the leakage jet the flow is in the middle of the vane passage and
in the next four vane passages, the flow is attached to the pressure surface of the passage and
separate from the suction surface.
Based on the above analysis, the guide vanes and the external turbine were modified
(second generation engine). The geometry is shown in Figure 65. The inlet blade angle of the
turbine was increased to be bigger than 90°, and the blade length was shortened to increase the
curvature. More guide vanes were added to narrow down the width of the guide vane passage to
avoid the flow separation and vortices. The outlet opening of the guide vane passage was made
smaller than the inlet opening to form a nozzle, so that the under-expanded exhaust gas can be
65

further expanded to the ambient pressure. The inlet and outlet of the turbine are expected to be at
the ambient pressure, and the kinetic energy change provides the power generation. The results
are listed in Table 4-20. It can be seen that the efficiency is 30.7% higher than the original
two-stage design.

Table 4-20: Results of the two-stage engine with a traditional external turbine
Cycle Power (kW) Efficiency (%)
2

8.2902

11.01

Figure 66 shows that stationary shock waves are generated when the outlet of the guide
vane passage and the inlet of the blade passage form a converging-diverging nozzle, which
causes a negative pressure gradient near the blade inlet generating positive torque.
Figure 67 shows that when an inner channel (channel A) opens to a guide vane passage
(passage B), the flow follows the guide vanes. However, when the channel wall (wall C) starts
passing the vane passage (passage D), the flow separates from the pressure surface of the passage;
when the inlet of a vane passage (passage E) is closed, the flow stops in the passage. Thus, the
flow is not continuous in the guide vane passages.

4.6.4.

Scalability

Section 4.5.1 states that the engine can be scaled without changing the ports timing and the
rotational speed should be scaled accordingly. In this section the second generation engine was

66

scaled down to 2/3 and the rotational speed was increased to 1.5 times. The results are listed in
Table 4-21. Comparing Table 4-21 and Table 4-20 it can be found that the efficiency maintains
and the power reduces to 2/3, which proves the scalability of the engine.
Table 4-21: Results of the scaled engine
Cycle Power (kW) Efficiency (%)
2

5.4813

10.92

67

CHAPTER 5. METHODS TO PREVENT BACKFIRE
For a WDE, not only does the leakage cause the energy loss, but also jeopardizes the safe
operation. Due to the pressure difference between the channels, the burned gas leaks from the
high pressure channels to the low pressure channels and the inlet ports. If the speed of the
leakage jet is too high, the jet will block the air from the air inlet and consequently destroy the
construction of the buffer layer between the burned gas and the fresh mixture and increase the
probability of the backfire. Besides, even though the gap is very small (0.3 mm) and the flame
cannot propagate through the gap, the leakage jet brings radicals and hot gas into the fuel inlet
and causes backfire. Thus, it is necessary to explore the method to avoid the leakage. This is
necessary when the premixed mixture is supplied from the fuel inlet. If only air goes through the
inlet port and the fuel is directly injected to the channels, backfire would not happen.

5.1. Simulation Methodology
Since the leakage flow is driven by the pressure difference, to provide a fare comparison the heat
addition was modeled by patch a high pressure and the according temperature to the channel
where the heat is released. In this analysis, 9 bar of the gauge pressure was used for the pressure
patching, and the patching temperature was calculated by the following equation

T3

9 ×105 Pa + 101325.1Pa ) T2
(
=
p2

(5.1)

The patching pressure is the same for all cases tested in this chapter and the working conditions
are the same in all cases in this chapter, listed in Table 5-1. The position of the patching channel
68

is shown in Figure 68. In the previous simulations, the channel is patched when its outlet is just
about to open so that the leakage would not develop much. In this analysis, the channel is
patched 1.5Tch before its outlet opens to give the leakage sufficient time to develop. However,
this patching method has some disadvantages. As the operation proceeds, the hot leakage gas will
concentrate on the outlet end of the suction surface. This gas has a very high temperature but is at
the same pressure as the rest area of the channel. After the pressure patching, the pressure rise
and the temperature rise are the same all over the channel, and the temperature in the hot leakage
area reaches a very high level. And then this hot gas leaks to the following channel. After the
next patching, its temperature gets higher again. As this process repeats, the temperature of the
hot leakage gas gets higher and higher. Its thermodynamic process is shown in Figure 69. For the
simulations, the temperature is limited under 5000 K, meaning that whenever the local
temperature gets above 5000 K, Fluent forces the temperature to be 5000 K while keeps other
variables unchanged. In this way, the temperature is reduced but the pressure does not change,
and consequently the density increases in this region, meaning that some additional mass is
added to this region. Fortunately, the temperature limiting occurs only near the outlet end of the
suction surface, so it does not affect the leakage in the inner gap. However, the efficiency is
affected.

69

Table 5-1: Working conditions of the simulations in CHAPTER 5

ptot,AI

0.1

(barg)
Ttot,AI (K)

ptot,FI

0.1

(barg)
300

Ttot,FI (K)

pGO

0

(barg)
300

TGO (K)

pTO

0

RPM

7500

(barg)
300

TTO (K)

300

Besides the pressure patching method, the temperature method was used and 2400 K was
used for the temperature patching. For the temperature patching method, no additional mass is
added. However, the temperature in the channel just before the patching is different in different
cases, so after the patching the pressure is different for different cases, which brings difficulty to
the comparison.
The third heat addition model is to specify a volumetric heat release rate to an area in the
channel by turning on the energy source term for a time period. The area is shown in Figure 70.
This area should not overlap with the hot leakage area near the end of the suction surface.
Otherwise due to the heat release in the hot leakage area, its temperature would go higher and
higher as in the pressure patching method. Moreover, the time for the heat release should be very
short. Since the heat release area is open, during the heat release the mass continues flowing out
of the area, causing a mass reduction in this area. Due to the constant heat release rate, the
specific heat released in this area increases with time and hence an increasing temperature rise
per time step. Eventually, the peak temperature will reach an unrealistic high value. In this
70

analysis, the heat release time was set to one time step. Thus, this method can be seen as an
energy patching method. If one channel is filled with the stoichiometric ethane-air mixture at 300
-5

K and 1 atm, the mass of ethane is 1.8194301×10

kg in the channel. The LHV of ethane is

6

47.80 MJ/kg, so the heat release is 47.80×10 ×1.8194301×10
channel is 2.6×10

-4

-5

= 869.7 J. The volume of the

3

m , so the volumetric heat release is 869.7/2.6×10
-7

the time step equal to 5×10

-4

3

= 3301909.2 J/m . For

s, the volumetric heat release rate is 3301909.2/5×10

-7

12

= 6.6×10

3

W/m . In this method, the cp was considered as a function of temperature. The values of cp at
different temperatures at 1 atm are tabulated in Table 5-2. A polynomial equation (Figure 71)
obtained from the six order curve fitting was used in Fluent to calculate cp.

Table 5-2: cp at different temperatures at 1 atm [19]

T (K)

cp (kJ/kg·K)

T (K)

cp (kJ/kg·K)

T (K)

cp(kJ/kg·K)

260

1005.3036

590

1048.8989

1400

1214.01

270

1005.3897

600

1051.1662

1450

1222.62

280

1005.5906

610

1053.4909

1500

1230.943

290

1005.8489

620

1055.8156

1550

1240.127

300

1006.1933

630

1058.169

1600

1249.024

310

1006.6238

640

1060.5511

1650

1258.495

320

1007.1117

650

1062.9332

1700

1267.966

330

1007.7144

660

1065.344

1750

1277.437

340

1008.3458

670

1067.7548

1800

1287.769

350

1009.092

680

1070.1656

1850

1298.388

71

Table 5-2 Con't
360

1009.8956

690

1072.5477

1900

1310.442

370

1010.7853

700

1074.9585

1950

1323.357

380

1011.8185

710

1077.3406

2000

1337.994

390

1012.8804

720

1079.7227

2050

1354.353

400

1014.0571

730

1082.1335

2100

1372.147

410

1015.2912

740

1084.5443

2150

1393.672

420

1016.6114

750

1086.9264

2200

1419.789

430

1018.0177

760

1089.2798

2250

1450.785

440

1019.5101

770

1091.6332

2300

1486.373

450

1021.0599

780

1093.9866

2350

1527.127

460

1022.6958

790

1096.34

2400

1573.908

470

1024.3891

800

1098.636

2450

1627.29

480

1026.1398

850

1110.403

2500

1688.134

490

1027.9479

900

1121.022

2550

1756.44

500

1029.8134

950

1131.928

2600

1836.8

510

1031.7363

1000

1141.973

2650

1925.77

520

1033.7453

1050

1151.731

2700

2026.22

530

1035.783

1100

1161.202

2750

2138.15

540

1037.8781

1150

1170.386

2800

2258.69

550

1040.0306

1200

1179.283

2850

2396.45

560

1042.1831

1250

1188.18

2900

2542.82

570

1044.393

1300

1197.077

2950

2697.8

580

1046.6316

1350

1205.687

3000

2858.52

To evaluate the effectiveness of each design, first a qualitative analysis will be performed by

72

examining the temperature contours and then a quantitative judgment will be done by comparing
the mass flow rates through the inlets with the no prevention design. If the mass flow rate is
higher than the no prevention design, the simulation using the energy patching method will be
performed to investigate the effect of the backfire prevention feature on the engine efficiency.

5.2. Various Backfire Prevention Designs
5.2.1.

No prevention design

To compare different backfire prevention designs, a no prevention design was tested without
any prevention feature. The geometry of the no prevention design is shown in Figure 68.
During the third operation cycle, the operation reaches the periodic mode. From the
temperature contour at the end of the third cycle shown in Figure 72, it can be seen that some hot
gas exists in the air inlet blocking the flow and the air layer fails to be created between the
exhaust gas and the mixture coming from the fuel inlet; some hot leakage gas exists at the inlet
corner of the channel that is about to close. The results from three patching method are listed in
Table 5-3.

5.2.2.

Plan A: U channel

Two U channels were added next to the fuel inlet and the air inlet respectively. The
configuration is shown in Figure 73. The idea of this feature is to redirect a part of the leakage jet
to counteract with the left of the jet. When the leakage jet passes the inlet of the U channel, a part
of is expected to flow through the channel. The jet out of the channel outlet and the jet passing

73

over the channel inlet are in different directions. When they merge, the direction of the merged
jet is expected to be along the channel walls instead of being in the circumferential direction. In
this way, the leakage flows into the rotor channels without causing a backfire in the fuel inlet or
blocking the inflow from the air inlet.
The temperature contour is shown in Figure 74. Compared with Figure 72, the leakage
prevention is not improved. The air inlet is still blocked and some hot gas still exists at the corner
of the fuel inlet. From Table 5-3, the mass flow rate decreases by 0.3398% in the pressure
patching simulation and increases by 1.812% in the temperature patching simulation. In either
case, the backfire prevention is ineffective.
From the velocity vector plot sequence of the left U channel during one channel period
(Figure ), the idea is realized when the channel wall is between the two arms of the U channel
and the left and the right arms are exposed to the high and lower pressure channels respectively.
When they are exposed to the same channel, the leakage jet neglects the U channel and passes
over it.
From the velocity vector plot sequence of the right U channel during one channel period
(Figure ), the same scenario happened to the right U channel.
Since this design is ineffective, the energy patching simulation was not performed to
compare the efficiency change.

5.2.3.

Plan B: modified U channel

If the distance between two arms of the U channel is greater than the channel inlet width,

74

the arms are always exposed to the different channels so that the inlet is always at a higher
pressure than the outlet. In this way the leakage gas can always flows through the U channel. The
U channel was modified based on this idea. The geometry of the modified U channel (MU
channel) is shown in Figure 77.
The vector plot sequence (Figure ) shows that the leakage gas flows through the left MU
channel continuously and the merged flow is in the direction of the outlet wall, which justifies
the idea of the MU channel. However, the backflow into the fuel inlet is observed as well.
The vector plot sequence (Figure ) shows that the leakage gas flows through the right MU
channel continuously and the merged flow is in the direction of the outlet wall. However, when
the jet out of the outlet impacts on the rotor channel wall, a part of it is deflected and goes into
the air inlet.
The temperature contour is shown in Figure 80. Compared with the previous designs, the
leakage to the air inlet does not change much, but that to the fuel inlet is more severe. From
Table 5-3, the mass flow rate decreases by 21.11% in the pressure patching simulation and
decreases by 10.92% in the temperature patching simulation. In either case, the leakage is worse.
From the pressure contour (Figure 81), it is found that the gauge pressure in the left MU
channel is between 0.855 bar and 0.765 bar which is much higher than that in the left U channel
which is between 0.270 bar and 0.225 bar. The gauge pressure in the right MU channel (0.135 ~
0.180 bar) is higher than that in the right U channel (0.045 ~ 0.090 bar). This is because the inlets
of the MU channels are connected to the regions at higher pressure. Although this makes the
leakage continuously flows through the channels, the high pressure merged flow causes more hot
75

gas to go into the inlets.
Since this design is ineffective, the energy patching simulation was not performed to
compare the efficiency change.

5.2.4.

Plan D: leakage pipe

In this design, the leakage jets on both sides of the inlets are taken away and re-injected to
the turbine to generate power. The configuration of a WDE with the leakage pipes is shown in
Figure 82.
The temperature contour is shown in Figure 83. It is observed that there is no hot gas in the
air inlet and a thin air layer is created behind the exhaust gas in the rotor channel; the hot gas
located at the corner of the fuel inlet and near the inlet of the rotor channel just closed is much
less than the no prevention design, which qualitatively prove the effectiveness of the leakage
pipes. From Table 5-3, the mass flow rate increases by 17.17% in the pressure patching
simulation and increases by 17.25% in the temperature patching simulation, which quantitatively
justifies this design. Then an energy patching simulation was performed to examine the effect of
the leakage pipes on the engine efficiency. However, even though the mass flow rate increases by
12.01%, the efficiency decreases by 1.88%.

5.2.5.

Plan E: jagged housing wall

In this design, some grooves are cut on the inner housing wall on both sides of the inlets to
bring more resistance to the leakage flow. The geometry is shown in Figure 84. The width of the
groove is 1.1 mm, the left wall is radial and the right wall has an angle of 30° with respect to the

76

periphery. There are 13 grooves on each side of the inlets.
From the temperature contour shown in Figure 85 it can be seen that the hot leakage
prevention is more effective than the leakage pipe design. From Table 5-3, the mass flow rate
increases by 19.25% in the pressure patching simulation and increases by 17.26% in the
temperature patching simulation. According to the result from the energy patching simulation,
the mass flow rate increases by 16.56% and the efficiency by 20.7%. Among all the backfire
prevention designs, the engine with the jagged housing wall has the highest mass flow rate and
the highest efficiency.

77

Table 5-3: Results of various backfire prevention designs
No

Leakage

Jagged

U

MU

prevention

pipe

wall

channel

channel

0.7357

0.8620

0.8773

0.7332

0.5804

0

17.17

19.25

-0.3398

-21.11

0.7671

0.8994

0.8995

0.7810

0.6833

0

17.25

17.26

1.812

-10.92

ppatching (barg)

7.759

8.246

8.891

7.802

6.652

 mean (kg/s)
m

0.7658

0.8578

0.8926

0

12.01

16.56

 mean (kg/s)
m

Pressure
patching

 mean
m

increment (%)
 mean (kg/s)
m

Temperature
patching

 mean
m

increment (%)

 mean
m

Energy

increment (%)

patching

power (kW)

27.47

26.93

33.19

Efficiency (%)

2.13

2.09

2.57

0

-1.88

20.7

Efficiency
increment (%)

78

-

CHAPTER 6. CONCLUSIONS AND FUTURE WORK
6.1. Summery and Conclusions
The principle of the WDE is described. Thermodynamic analysis of the WDE is performed on a
single channel and a full rotor. Through the single channel analysis, the general forms of the net
work extracted from the flow and the thermal efficiency are obtained in terms of the
mass-average properties. The results of the full rotor analysis are identical to those of the single
channel analysis. The analysis shows that under some certain conditions, the thermal efficiency
of the WDE can be equal to the Humphrey cycle efficiency. Thus, the mass-average flow can be
idealized as a Humphrey cycle.
A 1D numerical code is created using the TVDM scheme to predict the ports timing. The
wave diagram is obtained from the 1D code. The results from the 1D code and Fluent 2D
simulation show a good consistency between them.
The necessary geometric adjustments for a given rotor and working conditions are
introduced.
The energy loss of the WDE is examined by a 2D simulation of a single channel. The results
show that the two sources of loss are the significant outgoing kinetic energy and the
under-expansion of the burned gas.
The torque generation mechanism in the expansion process is investigated. For the WDE
channel, the torque is generated by the traveling expansion wave in the channel. The relationship

79

between the blade angle and the pressure gradient is obtained. Through the tests of various
channel shapes, it is discovered that smaller outlet opening of the channel leads to a longer
expansion process and hence a higher efficiency. The mechanism of the torque oscillation was
examined by the analysis of the pressure profiles. The analysis yields the conclusion that the
channel outlet should be at the outer radius. The convergent channel is considered as the best
channel shape and is implemented on a full WDE to improve the efficiency. With the convergent
channel the efficiency increases from 2.2% to 3.3%.
The effects of the peak pressure, the rotor dimension and the rotational speed on the engine
efficiency are investigated by the single channel simulations. The results show that a lower peak
pressure has a higher efficiency but a lower power. Whether to use a low peak pressure depends
on the priority of the power and efficiency. Bigger rotor dimension and higher rotational speed
can increase the efficiency by increasing the channel tip speed and consequently reducing the
outgoing velocity. The results are tested on the full engine for validation.
Two methods are proposed to reuse the energy from the exhaust gas – utilizing a return
channel and incorporating an external turbine. Preliminary investigation into the guide vanes and
the external turbine blades are performed. Incorporating an external turbine, the efficiency
increases to 11%.
To prevent the backfire, several plans are examined. By qualitative and quantitative analysis,
the best backfire prevention design is determined to be the jagged housing wall.

80

6.2. Suggestions for Future Work
The 1D code does not take into account any losses and 2D effect, such as the leakage, the heat
loss, the gradual opening and closing, the channel curvature and the variation of the channel
cross section area. In the future, the governing equation will be sophisticated to consider these
effects. In the current code, the outflow B.C. uses a simple extrapolation method. More advanced
outflow B.C. will be adopted.
In the channel shape investigation, only the expansion process is considered. However, the
shape best for the expansion may not be good for the scavenging and compression. In the future,
the channel performance for the full operation cycle will be studied. Moreover, more complex
shapes will be proposed to enhance the energy extraction performance.
Since the WDE with the external turbine has the best efficiency so far, more research will be
focused on the optimization of this design. The characteristics of the internal flow should be
investigated in detail. For example, due to the unsteady nature, the outlet condition of the
channel varies with time and consequently the inlet conditions of different nozzles are different
from each other. Since the nozzle passages are of the same shape, the nozzle outlet conditions are
different. Thus, as the blade passes the nozzles, its inlet sees different conditions. Therefore, the
blades operate on the off-design conditions for some time during one rotational cycle. In the
future, the nozzle passages should be redesigned to have different geometries such that the
nozzle outlets are at the identical condition.
The working range of the backfire prevention features should be examined. Combination of

81

various features may be more effective, which is of interest. Since the leakage is a very important
loss for wave engines, the research on the seal methods will continue, not only to avoid the
backfire but also to reduce the energy loss.
Since the efficiency of Humphrey cycle increases with the precompression ratio, methods to
enhance the precompression by the shock wave will be studied. Since the hammer shock wave is
very weak, a second shock wave can be generated by sending the burned gas back to the channel
with the mixture (Figure 86).
Besides, a critical factor is the flame speed. If the combustion takes a long time, the WDE
has to operate at low rotational speeds, which decreases the efficiency. More advance methods to
accelerate the combustion should be utilized. One of them is to use EGR to ignite the mixture.
In the future, with the high performance computers, 3D simulations can be carried out to
provide more realistic information.
This work concentrates on the numerical investigations only. Experiments should be done to
validate the conclusions of the numerical analysis.

82

APPENDIX

83

APPENDIX

Figure 1: Schematic model of a wave rotor inside the pressure wave supercharger (PWS) with its
development stages (soldered, cast with single of channels, cast with double row of channels)
and actual NASA four port wave rotor drum (right). For interpretation of the references to color
in this and all other figures, the reader is referred to the electronic version of this dissertation.

84

Figure 2: Porting and flow scheme for a through-flow (left) and reverse flow wave rotor (taken
from [55])

Figure 3: Assembly of successfully working wave engine with wave rotor (left) and rear and
front stator plates [3]

85

Figure 4: ABB test rig (left) and wave rotor combustor cross section (right) [4]

expansion
work
Figure 5: By using shockwaves for internal work transmission, the Wave Disk Engine virtually
eliminates the hardware (shown in gray) in gas turbines and internal combustion engines that
serve to transmit expansion work internally to perform the needed inlet gas compression.

86

1/5

4

2

3

Figure 6: Schematic of an internal combustion WDE

Figure 7: T-s diagram of the mass-average flow in a WDE

87

Reflected
compression
wave

Exhaust gas leaving

Hammer
shock wave

x

t

Incident
expansion wave

Fresh mixture entering

Reflected
expansion wave
(a) Pressure contour

Interface between the exhaust
gas and the fresh mixture

Exhaust gas leaving

x

t

Fresh mixture entering
(b) Temperature contour
Figure 8: Wave diagram obtained from the TVDM scheme

88

5

7

7

x 10

1.2
1.2

Pressure (bar)

1.1

6

Pressure(Pa)

Density(kg/m3)

6

55

11

3

Density (kg/m )

0.9

0.8
0.8

44

0.7

0.6
0.6

33

0.5

500

400
400
300

200
200

Temperature(K)

Velocity(m/s)

0.4
0.4
2
20.05
0.05
0.06
0.06
0.07
0.07
0.08
0.08
0.09
0.09
0.1
0.1
0.05 0.06
0.07
0.05
0.06
0.07
x(m)0.08 0.09 0.1
x(m) 0.08 0.09 0.1
x (m)
x (m)
2100
800
800
2100
Temperature (K)
700
2000
Velocity (m/s)
600
600
1900
1900
1800

1700
1700
1600

1500
1500
100
1400
00
1400
0.05
0.05
0.06
0.06
0.07
0.08
0.08
0.09
0.09
0.1
0.1
0.05
0.06
0.07
0.05
0.06 0.07
0.07
x(m) 0.08 0.09 0.1
x(m) 0.08 0.09 0.1
x (m)
x (m)

Figure 9: Profiles of pressure, density, velocity and temperature for TVDM scheme (black) and
GT scheme (blue) at 5.44e-5 s and its corresponding location in the wave diagram

89

4.5
4

0.8

Pressure (bar)

3

Density (kg/m )
0.7

3.5
0.6

3
2.5

0.5
2
0.45
0.05 0.06 0.07 0.08 0.09 0.1
0.05 0.06 0.07 0.08 0.09 0.1
x (m)
x (m)
800
1800
Velocity (m/s)
Temperature (K)
600
1700
400

1600

200

1500
1400
0
0.05 0.06 0.07 0.08 0.09 0.1
0.05 0.06 0.07 0.08 0.09 0.1
x (m)
x (m)

Figure 10: The profiles of pressure, density, velocity and temperature from GT scheme (blue) and
TVDM scheme (black) at 7e-5 s and its corresponding location in the wave diagram

90

1.1

1.2

Pressure (bar)

1

1

0.9

0.8

0.8

0.6

0.7

0.4

3

Density (kg/m )

0.6
0.2
0.05 0.06 0.07 0.08 0.09 0.1
0.05 0.06 0.07 0.08 0.09 0.1
x (m)
x (m)
300
1200
Velocity (m/s)
250
1000
200
800
150

600

100

400

50
0.05 0.06 0.07 0.08 0.09 0.1
x (m)

Temperature (K)

200
0.05 0.06 0.07 0.08 0.09 0.1
x (m)

Figure 11: Profiles of pressure, density, velocity and temperature from (b) GT scheme (blue) and
TVDM scheme (black) at 2.3e-4 s and its corresponding location in the wave diagram (a)

91

2

2
1.8

1.6

1.6
1.4

1.2

1.2

Pressure (bar)

0.8
0.05 0.06 0.07 0.08 0.09 0.1
x (m)
190
Velocity (m/s)
150

3

Density (kg/m )

1
0.05 0.06 0.07 0.08 0.09 0.1
x (m)
360
340

110

320

70

300

30

Temperature (K)
280
0.05 0.06 0.07 0.08 0.09 0.1
x (m)

-10
0.05 0.06 0.07 0.08 0.09 0.1
x (m)

Figure 12: Profiles of pressure, density, velocity and temperature from (b) GT scheme (blue) and
TVDM scheme (black) at 5.33e-4 s and its corresponding location in the wave diagram (a)

92

3

Specific Volume (kg/m )
Figure 13: T-s diagram (upper) and p-v diagram (lower) for one fluid particle

93

α_EO2

EO_closed
α_AI1
α_AI2

α_FI1
α_FI2

α_EO1
EO_open
FI_closed
FI_open
AI_closed
AI_open

Gap

Inner diameter
β1

β4

Outer diameter

Figure 14: Geometric parameters of a WDE

94

Figure 15: Heater for ignition

8.74e4

2600

0

281
Pressure contour (Pa)

Temperature contour (K)

Figure 16: Backflow into the air inlet

95

0.055

0
(b)

(a)

Figure 17: Relative velocity vector plot (a) and contour of mass fraction of methane (b)

(a)

(b)

Figure 18: Absolute velocity vector plot near the inlet (a) and outlet (b)

96

3440

0.055

280

0
Temperature contour (K)

Contour of mass fraction of methane

Figure 19: Backfire into the fuel inlet

97

0.055

0

Contour of mass fraction of methane
Figure 20: Incomplete combustion

98

2000

0.055

287

0
(a)

(b)

(c)
Figure 21: (a) temperature contour, (b) mass fraction of methane contour and (c) relative velocity
vector plot for the baseline design

99

Figure 22: Power history in the second cycle

100

T

s
Figure 23: T-s diagrams of 10 particles in the channel

101

Outlet posts

Fuel inlet posts
Air inlet posts

Channel

Specific enthalpy

Pressure
p (Pa)

Figure 24: Geometry for the 2D investigation

Time

Time

Figure 25: Histories of the outlet enthalpy (left) and outlet pressure (right) of the 2D simulation

102

Figure 26: Geometry of the C2 channel

8.55e-3

-1.63e+3
Figure 27: Pressure contours (Pa) in the curved channels

103

8
6

Expansion

Compression

4
Cm

Scavenging
2
0
-2
-4
0.015

0.0155

0.016

0.0165
t (s)

0.017

0.0175

Figure 28: History of the moment coefficient on the channel walls

Figure 29: Geometry of the baseline shape

104

0.018

M0

Figure 30: Profiles of a shock wave (top) and a planar acoustic wave (bottom) through a
two-dimensional right-angle bend, numbers 1 - 6 are representative wave front configurations
within the bend (Edwards et al [12], Kentfield [23])

105

6.1e5

5.5e5
(1)

(2)

(3)

(4)

(6)
(5)
Figure 31: Sequence of the pressure contours when the channel just opens

106

6.1e5

3.0e5
Figure 32: Sequence of the pressure contours during the incident expansion wave propagates

4.74e5

3.69e5

1.00e5

1.00e4
(1)

(2)

2.93e5

2.34e5

-3.0e4

-5.0e4
(3)

(4)

Figure 33: Sequence of the pressure contours when the reflecting wave interferes with the
incident wave

107

6e5
a
b
0
Figure 34: Pressure contours as the expansion wave propagates in a straight channel

Figure 35: Pressure contour of the C2 channel

Figure 36: Geometry of the C3 channel

108

Figure 37: Geometries of J1 (left) and J2 (right) channels

Efficiency (%)

4
3.5
3
2.5
2
0.25

0.3

0.35
0.4
0.45
Expansion Duration (ms)

0.5

Figure 38: Efficiency vs. expansion duration for C1, C2, C3, J1 and J2 channels

Figure 39: Outlet opening (circle radius)

109

Expansion Duration (ms)

0.5
0.45
0.4
0.35
0.3
0.25
0.4

0.5

0.6
0.7
Outlet Opening (cm)

0.8

Figure 40: Expansion duration vs. outlet opening

Figure 41: Geometries of the MJ2(left) and MC3(right) channels

110

Figure 42: Geometry of the convergent channel

111

8
7
6
C1

Cm

5

C2

4

C3

3

J1

2

J2

1
0
-1
0.00E+00

1.00E-04

2.00E-04

3.00E-04

4.00E-04

5.00E-04

t (s)
8
7
6

Cm

5
4

MJ2
MC3
convergent

3
2
1
0
-1
0.00E+00

4.00E-04

8.00E-04
t (s)

1.20E-03

1.60E-03

Figure 43: Torque histories of five short expansion cases (upper) and three long expansion cases
(lower)

112

7e5
6e5
5e5
4e5
Static
3e5
Pressure
2e5
(Pa)
1e5
0
-1e5

Positive torque

5
7e5
6e5
5e5
4e5
Static
Pressure 3e5
2e5
(Pa)
1e5
0
-1e5

6
7
8
9
10
Radial Coordinate (cm)

Positive torque

5
7e5
6e5
5e5
Static
4e5
Pressure 3e5
(Pa)
2e5
1e5
0
-1e5

6
7
8
9
10
Radial Coordinate (cm)

Positive torque

5

6
7
8
9
10
Radial Coordinate (cm)

Figure 44: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the propagation of the incident
expansion wave

113

7e5
6e5
5e5
4e5
Static
Pressure 3e5
2e5
(Pa)
1e5
0
-1e5

Positive torque

5
Figure 44 Con't

114

6
7
8
9
10
Radial Coordinate (cm)

5e5

0
7e5
6e5
5e5
4e5
3e5
2e5
1e5
0
0.09

0.092

0.094

0.096

0.098

0.1

Figure 45: Pressure contour (top) and pressure profiles (bottom) near the trailing edge of the
channel walls when the channel outlet is not fully open

115

5e5

-1e5
6e5
5e5
4e5
3e5
2e5
1e5
0
-1e5
0.09

0.092

0.094

0.096

0.098

0.1

Figure 46: Pressure contour (top), pressure profiles (middle) and Mach number contour near the
trailing edge when the channel outlet is fully open

116

3

1
Figure 46 Con't

117

6e5
5e5
4e5
3e5
Static
Pressure 2e5
1e5
(Pa)
0
-1e5
-2e5

Negative torque

Positive torque

5
6e5
5e5
4e5
Static
3e5
Pressure 2e5
(Pa)
1e5
0
-1e5
-2e5

6
7
8
9
10
Radial Coordinate (cm)

Negative torque
Positive torque

5
6e5
5e5
4e5
3e5
Static
2e5
Pressure
1e5
(Pa)
0
-1e5
-2e5

6
7
8
9
10
Radial Coordinate (cm)

Negative torque

Positive torque

5

6
7
8
9
10
Radial Coordinate (cm)

Figure 47: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the propagation of the reflected
expansion wave

118

Figure 48: Relative velocity vectors showing the back flow near the inlet end

119

5e5
4e5
3e5
2e5
Static
Pressure 1e5
0
(Pa)
-1e5
-2e5
5

6
7
8
9
10
Radial Coordinate (cm)

5

6
7
8
9
10
Radial Coordinate (cm)

5e5
4e5
3e5
Static
2e5
Pressure
1e5
(Pa)
0
-1e5
-2e5

Figure 49: Sequences of the pressure contour (left )and the pressure distribution (right) on the
pressure surface (black) and suction surface (red) showing the positive pressure gradient caused
by the back flow near the inlet

120

Figure 50: Pressure contour in the beginning of the simulation

121

cycle 1
100
80

Cm

60
40
20
0
-20
-40
0.00E+00

5.00E-04

1.00E-03

1.50E-03
t (s)

2.00E-03

2.50E-03

3.00E-03

0.005

0.0055

0.006

cycle 2
100
80

Cm

60
40
20
0
-20
-40
0.003

0.0035

0.004

0.0045
t (s)

Figure 51: Torque histories of the first three cycles

122

cycle 3
100
80

Cm

60
40
20
0
-20
-40
0.006

0.0065

0.007

0.0075
t (s)

0.008

0.0085

Figure 51 Con't

Figure 52: Geometry of the single rotor engine with convergent channels

123

0.009

Cm

cycle 1
70
60
50
40
30
20
10
0
-10
-20
0.00E+00

5.00E-04

1.00E-03

1.50E-03
t (s)

2.00E-03

2.50E-03

3.00E-03

0.005

0.0055

0.006

Cm

cycle 2
70
60
50
40
30
20
10
0
-10
-20
0.003

0.0035

0.004

0.0045
t (s)

Figure 53: Histories of torque of the first two cycles

124

5

10

p = 2 bar >1bar
Flow is under-expanded.

Figure 54: Profile of pressure in case 1 when the incident expansion propagates

5

10

p = 1bar
Flow is fully expanded.

Figure 55: Profile of pressure in case 2 when the incident expansion propagates

125

Figure 56: WDE with C3 channels

wr = cr
c
cθ

wθ
U

w

Figure 57: Velocity triangle at the trailing edge of a channel wall

126

11

Efficiency (%)

10
9
8
7
6
5
4
3
10

15
Outer Radius (cm)

20

Efficiency (%)

Figure 58: Efficiency vs. outer radius

10
9
8
7
6
5
4
3
2
1
0
1000

5000

10000
RPM

15000

Figure 59: Efficiency vs. RPM

127

20000

Figure 60: Geometry of the engine with 30 cm OD

128

Figure 61: WDE with a return channel

129

Figure 62: Configuration of a two-stage engine

130

5e4

B
A

-1e4
Figure 63: Pressure contour in the external turbine

131

Figure 64: Streamlines in the guide vane passage

132

Figure 65: Configuration of a second generation engine

133

5e4

-6e4
Figure 66: Pressure contour in the blade passages of the second generation engine

134

C
A
D
B

E
Figure 67: Streamlines in the guide vane passage of the second generation engine

135

Figure 68: Temperature contour of the no prevention design at 2e-5 s

136

Pmax

Po

Figure 69: Thermodynamic process of the hot leakage gas

137

Figure 70: Energy patching area (red area)

138

3000

y = 6e-14x5 - 2e-10x4 + 3e-08x3 + 0.0003x2 - 0.081x + 1002.1

2800
2600

Experimental data
Fiitted curve

cp (kJ/kg·K)

2400
2200
2000
1800
1600
1400
1200
1000
0

500

1000

1500
T (s)

2000

2500

3000

Figure 71: Variation of cp with T of the experimental data and the fitted curve

139

1000

290
Figure 72: Temperature contour at the end of the third cycle from the pressure patching method
for the no prevention design

140

Figure 73: Idea of the U channel

141

1000

290
Figure 74: Temperature contour at the end of the third cycle from the pressure patching method
for the U channel design

(1)
(2)
Figure 75: Velocity vector plot sequence of the left U channel during one channel period

142

(3)

(4)

(5)

(6)

(7)

(8)

Figure 75 Con't

143

(1)

(2)

(3)
(4)
Figure 76: Velocity vector plot sequence of the right U channel during one channel period

144

(5)

(6)

(7)

(8)

Figure 76 Con't

145

Figure 77: Geometry of the MU channel

(1)
(2)
Figure 78: Velocity vector plot sequence of the left MU channel during one channel period

146

(3)

(4)

(5)

(6)

(7)

(8)

Figure 78 Con't

147

(1)

(2)

(3)
(4)
Figure 79: Velocity vector plot sequence of the right MU channel during one channel period

148

(5)

(6)

(7)

(8)

Figure 79 Con't

149

1000

290
Figure 80: Temperature contour at the end of the third cycle from the pressure patching method
for the MU channel design

150

9000

0

9000

0
Figure 81: Comparison between the pressure in the MU channels (top) and U channels (bottom)

151

Figure 82: Configuration of a WDE with the leakage pipes

152

1000

290
Figure 83: Temperature contour at the end of the third cycle from the pressure patching method
for the leakage pipe design

153

Figure 84: Geometry of the jagged housing wall

154

1000

290
Figure 85: Temperature contour at the end of the third cycle from the pressure patching method
for the jagged wall design

155

0.055

0

Contour of mass fraction of methane
Figure 86: Using EGR to generate a shock wave and ignition

156

REFERENCES

157

REFERENCES
[1] Anderson, John, 1995, “Computational Fluid Dynamics: the Basics with Applications”,
McGraw-Hill Science/Engineering/Math
[2] Aichlmayr, H. T., Kittelson, D. B., and Zachariah, M. R., “Micro-HCCI Combustion:
Experimental Characterization and Development of a Detailed Chemical Kinetic Model with
Coupled Piston Motion,” Combustion and Flame, Vol. 135, No. 3, 2003, pp. 227-248.
[3] Akbari, P., Nalim, M. R., Mueller, N., “A Review of Wave Rotor Technology and Its
Applications,” ASME Journal of Engineering for Gas Turbines and Power, October 2006,
Vol. 128, pp. 717-735.
[4] Akbari, P. and Nalim, M. R., “Review of Recent Developments in Wave Rotor Combustion
Technology.” Journal of Propulsion and Power, Vol. 25, No. 4, pp 833-844, July–August
2009.
[5] Akbari, P., and Nalim, M. R., “Analysis of Flow Processes in Detonative Wave Rotors and
Pulse Detonation Engines,” AIAA Paper 2006-1236, 2006.
[6] Akbari, P., Alparslan, B., Kilchyk, V., and Nalim, M. R., “Numerical Analysis of
Hydrogen-Fueled Wave Rotors for Gas Turbine Applications,” First International Hydrogen
Energy Congress, Turkey, 2005.
[7] Azim A., 1974, “An Investigation Into the Performance and Design of Pressure Exchangers,”
Ph.D. Thesis, University of London.
[8] Azoury P.H., Hai S.M., 1975, Computerized analysis of dynamic pressure exchanger
scavenge process, Proc. I Mech E vol. 189/18/75 pp. 149-158
[9] Board on Army Science and Technology, “Meeting the Energy Needs of Future Warriors”
National Research Council of The National Academies, Washington, DC, 2004.
[10] Brouillette, M., “ Shock Waves at Microscales,” Shock Waves, Vol. 13, No. 1, 2003, pp.
3-12.
[11] Clarke, J. F., “Fast Flames, Waves and Detonation,” Progress in Energy and Combustion
Science, Vol. 15, No. 3,1989, pp. 241-271.
[12] Edwards, D. H., Fearnley, P. and Nettleton, M. A., “Shock diffraction in channels with 90
deg bends”, Journal of Fluid Mechanics (ISSN 0022-1120), vol. 132, July 1983, p. 257-270.

158

[13] Epstein, A. H., “Millimeter-Scale, Micro-Electro-Mechanical Systems Gas Turbine Engines,”
ASME Journal of Engineering for Gas Turbines and Power, Vol. 126, No. 2, 2005, pp.
205-226.
[14] Fatis A., and Ribaud Y., 1997, “Numerical Analysis of the Unsteady Flow Inside Wave
Rotors Applied to Air Breathing Engines,” 13 International Symposium on Airbreathing
Engines
[15] Fatis A., Lafond and Ribaud Y., 1998, “Preliminary Analysis of the Flow Inside a
Three-Port Wave Rotor by Means of a Numerical Model,” Aerospace Science and
Technology, Vol. 2, No. 5, pp. 289-300
[16] Fatis A., and Ribaud Y., 1999, „Thermodynamic Analysis of Gas Turbines Topped with
Wave Rotors,” Aerospace Science and Technology, Vol. 3, No. 5, pp. 293-299.
[17] Fernandez-Pello, A. C., Pisano, A. P., Fu , K., Walther, D. C., Knobloch, A. K., Martinez, F.
C., Senesky, M., Stoldt, C., Maboudian, R. M., Sanders, S., and Liepmann, D., “MEMS
Rotary Engine Power System,” Proceedings of IEEE Transactions on Sensors and
Micromachines, Vol. 123, No. 9, 2003, pp. 326-330.
[18] Fong, K. K., and Nalim, M. R., “Gas Dynamic Limits and Optimization of Pulsed
Detonation Static Thrust,” AIAA Paper 2000-3471, 2000.
[19] Hilsenrath, Joseph; et al.: Tables of Thermal Properties of Gases. Circ. 564, National Bureau
of Standards, Nov. 1, 1955
[20] Iancu. F, Piechna. J, Dempsey. E, and Müller. N, “The Ultra-Micro Wave Rotor Research a t
Michigan State University, second International Symposium on Innovative Aerial/Space
Flyer Systems, The University of Tokyo, 2005, pp 65–70.
[21] Iancu, F., Piechna. J., Dempsey. E and Müller N, “Ultra-Micro Wave Rotor Investigations,”
fifth International Workshop on Micro and Nanotechnology for Power Generation and
Energy Conversion Applications, Japan, 2005, pp 93–96.
[22] Izzy, Z., and Nalim, M. R., “Wave Fan and Rotary-Ejector Pulsed Performance Prediction,”
AIAA Paper 2002-4068, 2002.
[23] Jonsson, V. K., Matthews, L., and Spalding, D. B., 1973, “Numerical Solution Procedure for
Calculating the Unsteady, One-Dimensional Flow of Compressible Fluid,” ASME Paper
73-FE-30.
[24] Kentfield, J. A. C., 1993, Nonsteady, one-dimensional, internal, compressible flows, Oxford
University Press

159

[25] Kentfield, J. A. C, Thermodynamics of airbreathing pulse-detonation engines, Journal of
Propulsion and Power (2002), Volume: 18, Issue: 6, Publisher: American Inst. Aeronautics
and Astronautics Inc., Pages: 1170-1175
[26] Litke, P. J., Schauer, F. R., Paxson, D. E., Bradley, R. P., and Hoke, J. L., “Assessment of
the Performance of a Pulsejet and Comparison with a Pulsed-Detonation Engine,” AIAA
Paper 2005-0228, 2005.
[27] Matthews L., 1969, “An Algorithm for Unsteady Compressible One-Dimensional Fluid
Flow,” M.S. Thesis, University of London.
[28] Mawid, M. A., and Sekar, B., “A Numerical Study of Active Control of Combustion-Driven
Dynamic Instabilities in Gas-Turbine Combustors,” AIAA Paper 1999-2778, 1999.
[29] Nalim, M. R., “Longitudinally Stratified Combustion in Wave Rotors,” Journal of
Propulsion and Power, Vol. 16, No. 6, 2000, pp. 1060-1068.
[30] Nalim, M. R., and Paxson, D. E., “A Numerical Investigation of Premixed Combustion in
Wave Rotors,” Journal of Engineering for Gas Turbines and Power, Vol. 119, No. 3, 1997,
pp. 668-675.
[31] Nalim, M. R., and Jules, K., “Pulse Combustion and Wave Rotors for High-Speed
Propulsion Engines,” AIAA Paper 98-1614, 1998.
[32] Nalim, M. R., and Izzy, Z, “Rotary Ejector Enhanced Pulsed Detonation System,” AIAA
Paper 2001-3613, 2001.
[33] Nalim, M. R., and Izzy, Z, “Simulation of a Wave Rotor Pulse Detonation Engine with
Integrated Ejector,” ISABE Paper 2001-1214, 2001.
[34] Nalim, M. R., “Rotary Ejector Enhanced Pulsed Detonation System and Method,” US Patent
6,845,620, 2005.
[35] Paxson D. E.,1992, "A General Numerical Model for Wave Rotor Analysis,"NASA TM
105740, July.
[36] Paxson D. E., 1993, „A Comparison Between Numerically Modeled and Experimentally
Measured Loss Mechanisms in Wave Rotors,” AIAA Paper 93-2522.
[37] Paxson D. E., 1995,"A Comparison Between Numerically Modelled and Experimentally
Measured Loss Mechanisms in Wave Rotors," AIAA Journal of Propulsion and Power, Vol.
11, No. 5, pp. 908-914, (also NASA TM 106279).

160

[38] Paxson D. E., 1995, “A Numerical Model for Dynamic Wave Rotor Analysis,” AIAA Paper
95-2800. Also NASA TM-106997.
[39] Paxson D. E.,1996, "Numerical Simulation of Dynamic Wave Rotor performance," AIAA
Journal of Propulsion and Power, Vol. 12, No. 5, pp.949-957, (also, NASA TM 106997).
[40] Paxson D. E., 1997, "A Numerical Investigation of the Startup Transient in a Wave Rotor,"
ASME Journal of Engineering for Gas Turbines and Power, Vol.119, No. 3, pp. 676-682,
also ASME Paper 96-GT-115, June, 1996, (also, NASA TM 107196).
[41] Paxson D. E., 1997, “A Numerical Investigation of the Startup Transient in a Wave Rotor,”
Journal of Engineering for Gas Turbines and Power, 119, No. 3, pp. 676-682. Also ASME
Paper 96-GT-115, and NASA TM 107196.
[42] Paxson, D. E., 1998, “An Incidence Loss Model for Wave Rotors With Axially Aligned
Passages, ” AIAA-98-3251, July, 1998; also NASA TM-1998-207923.
[43] Paxson, D. E., “Wave Augmented Diffusers for Centrifugal Compressors,” AIAA Paper
98-3401, 1998.
[44] Paxson, D. E., “An Improved Numerical Model for Wave Rotor Design and Analysis,”
AIAA Paper 93-0482, 1993.
[45] Paxson, D. E., and Wilson, J., “Recent Improvements to and Validation of the One
Dimensional NASA Wave Rotor Model,” NASA TM-106913, 1995.
[46] Paxson, D. E., “A Numerical Investigation of the Startup Transient in a Wave Rotor,”
Journal of Engineering for Gas Turbines and Power, Vol. 119, No. 3, 1997, pp. 676-682.
[47] Paxson, D. E., and Lindau, J. W., “Numerical Assessment of Four-Port Through-Flow Wave
Rotor Cycles with Passage Height Variation,” AIAA Paper 97-3142, 1997.
[48] Paxson, D. E., “A Sectored-One-Dimensional Model for Simulating Combustion
Instabilities in Premix Combustors,” AIAA Paper 2000-0313, 2000.
[49] Paxson, D. E., and Perkins, H. D., “Thermal Load Considerations for Detonative
Combustion-Based Gas Turbine Engines,” AIAA Paper 2004-3396, 2004.
[50] Paxson, D. E., “Performance Evaluation Method for Ideal Airbreathing Pulse Detonation
Engines,” Journal of Propulsion and Power, Vol. 20, No. 5, 2004, pp. 945-947.
[51] Perkins, H. D., Paxson, D. E., Povinelli, L. A., Petters, D. P., Thomas, S. R., Fittje, J. E., and
Dyer, R. S., “An Assessment of Pulse Detonation Engine Performance Estimation Methods
Based on Experimental Results,” AIAA Paper 2005-3831, 2005.
161

[52] Piechna J 2005 Wave machines, models and numerical simulation Oficyna Wydawnicza
Politechniki Warszawskiej,Warszawa
[53] Piechna, J. R., “Feasibility Study of the Wave Disk Micro-Engine Operation,” Journal of
Micromechanics and Microengineering, Volume 16, Number 9, pp. 270-281, 2006.
[54] Pohořelský, L., Obernesser, P., Vavra J., Klir, V., Macek, J, “1-D Model and Experimental
Tests of Pressure Wave Supercharger.” ASME Paper IMECE2007-43427, Seattle 2007,
Washington, USA.
[55] Pohořelský, L., Sané P. A., Rozsas T., Müller, N, "Wave Rotor Design Procedure for Gas
Turbine Engine Enhancement." Proceedings of ASME Turbo Expo 2008, Berlin, June 2008.
[56] Resler E., Moscari J., Nalim R., 1994, Analytic design methods for wave rotor cycles, AIAA
Joumal of Propulsion and Power 10 (5)
[57] Spadaccini, C. M., Mehra, A., Lee, J., Zhang, X., Lukachko, S., and Waitz, I. A., “High
Power Density Silicon Combustion Systems for Micro Gas Turbine Engines,” ASME Journal
of Engineering for Gas Turbines and Power, Vol. 125, No. 3, 2003, pp. 709-719.
[58] Spadaccini, Peck, J., and Waitz, I. A., “Catalytic Combustion Systems for Microscale Gas
Turbine Engines,” ASME Journal of Engineering for Gas Turbines and Power, Vol. 129, No.
1, 2007, pp. 49-60.
[59] Snyder, P., Alparslan, B., and Nalim, M. R., “Gas Dynamic Analysis of the Constant
Volume Combustor, A Novel Detonation Cycle,” AIAA Paper 2002-4069, 2002.
[60] Snyder, P., Alparslan, B., and Nalim, M. R., “Wave Rotor Combustor Test Rig Preliminary
Design,” ASME Paper IMECE2004-61795, 2004.
[61] Horrocks G.D., Reizes J.A., Mallinson S.G., 1998, Numerical study of flow in a rotating
shock tube, 13th Australasian Fluid Mechanics Conference Monash University, Melbourne,
Australia 13-18 December 1998
[62] MacCormack, Robert W.: The Effect of Viscosity in Hypervelocity Impact Cratering.
AIAA-69-354, Apr.-May 1969.
[63] D Gottlieb and E Turkel, Dissipative two–four methods for time dependent problems. Math.
Comp., 30 (1976), p. 703
[64] Dongfang Liang, BinLiang Lin, Roger Alexander Falconer, Simulation of rapidly varying
flow using an efficient TVD-MacCormack scheme, International Journal for Numerical
Methods in Fluids (2007) Volume: 53, Issue: August 2006, Pages: 811-826

162

[65] William H Heiser, David T Pratt, Thermodynamic cycle analysis of pulse detonation engines,
Journal of Propulsion and Power (2002), Volume: 18, Issue: 1, Publisher: American Inst.
Aeronautics and Astronautics Inc., Pages: 68-76
[66] E. Wintenberger and J. E. Shepherd "Thermodynamic Analysis of Combustion Processes for
Propulsion Systems." 4second AIAA Aerospace Sciences Meeting and Exhibit, Jan 5-8,
Reno NV 2004. AIAA2004-1033.
[67] Weber, Helmut E., 1995, Shock Wave Engine Design, John Wiley and Sons, New York
[68] Wilson J. and Paxson D.E.: Jet engine Performance Enhancement through Use of a
Wave-Rotor Topping Cycle. NASSA TM-4486, 1993
[69] Zehnder G., 1970,Berechnungsaufgaben bei der Entwicklung des Comprex.Schweizerische
Bauzeitung, 88. Jahrgang, Heft 30, 23. Juli, 1970.
[70] Zehnder G., 1971, Berechnung von Druckwellen in der Auf ladetechnik.Brown Boveri Mitt
4/5 - 71.
[71] Zehnder G., 1971, “Calculating Gas Flow in Pressure-Wave Machines”Brown Boveri
Review, No. 4/5, pp. 172-176.

163