Cranfield University

Brahmadevan V Padmarajan

Numerical Modelling and Simulation of Rotary Engine

School of Engineering

MSc

Cranfield University
School of Engineering

MSc Thesis

2004

Brahmadevan V Padmarajan

Numerical Modelling and Simulation of Rotary Engine

Supervisor: Dr. Philip Rubini

Academic Year 2003 to 2004

This thesis is submitted in partial fulfilment of the requirements

for the degree of MSc
Cranfield University, 2004. All rights reserved. No part of this publication may be
reproduced without the written permission of the copyright holder.

ABSTRACT
Smooth, powerful, compact and light are the salient features of rotary engine. For many
peoples surprise rotary engines are still used in many field ranging from military to
electronics not to forget mazda RX-8.
In order to design efficient rotary engines, it is necessary to have a good understanding
of how engine design and operating parameters affect the unsteady multidimensional
fluid flow, fuel-air mixing and combustion occurring inside rotary engine.
Mathematical models have been developed in the last couple of decades to enhance
understanding of the flow field, fuel-air mixing and combustion in rotary engine. Earlier
improvements in performance and emission have been achieved primarily through
experimental trial and error procedures making the analysis time consuming and costly.
FLUENT 6.1 was used to develop a flexible rotary engine model. Dynamic mesh
capabilities are explored in handling the eccentric rotation of the rotor and deforming
housing of the rotary engine.
Simulation is carried out with the model having no inlet and exhaust ports and rotor
pocket for demonstration and finally recommendation and suggestions are made for
future work.

ACKNOWLEDGMENT
At the end of my thesis I would like to thank all those people who made this thesis
possible and an enjoyable experience for me.
I am deeply indebted to my supervisor, Dr. Philip Rubini, for patient guidance,
encouragement and excellent advice. Without his help, this work would not be possible.
I am grateful to Dr. Jin Yan for his generous assistance during this time.
Finally, thanks to my parents and numerous friends for always offering support and
love.
I dedicate this thesis to my father and mother.

ii

TABLE OF CONTENTS
ABSTRACT ...................................................................................................................... i
ACKNOWLEDGMENT .................................................................................................. ii
TABLE OF FIGURES ..................................................................................................... v
TABLE OF EQUATIONS ............................................................................................. vii
TABLE OF SYMBOLS AND NOTATIONS............................................................... viii
1
Introduction .............................................................................................................. 1
2
Review of Rotary engine .......................................................................................... 3
2.1
History of Rotary engine .................................................................................. 3
2.2
Principle of Operation ...................................................................................... 3
2.3
Salient features of rotary engine....................................................................... 4
2.4
Basic Equation of the rotary engine ................................................................. 5
2.5
Basic construction of the rotary engine ............................................................ 7
2.5.1
Rotor ......................................................................................................... 7
2.5.2
Housing..................................................................................................... 7
2.5.3
Phasing gear mechanism .......................................................................... 8
2.5.4
Output shaft .............................................................................................. 8
2.5.5
Gas sealing mechanism ............................................................................ 8
2.5.6
Intake and exhaust port system................................................................. 9
2.6
RENESIS .......................................................................................................... 9
2.7
Other applications of rotary engine ................................................................ 10
2.8
Scope for improvement .................................................................................. 12
2.9
Conclusion ...................................................................................................... 12
3
Review of the analysis of rotary engine using computational fluid dynamics....... 13
3.1
Computational Fluid Dynamics...................................................................... 13
3.2
Review of the analysis of RE using CFD ....................................................... 13
3.3
Conclusion ...................................................................................................... 17
4
Dynamic mesh ........................................................................................................ 19
4.1
Introduction to dynamic mesh ........................................................................ 19
4.2
Dynamic mesh update methods...................................................................... 19
4.2.1
Spring based smoothing ......................................................................... 19
4.2.2
Dynamic layering ................................................................................... 21
4.2.3
Local remeshing method ........................................................................ 23
4.3
Volume mesh update procedure ..................................................................... 25
4.4
Solid-Body Kinematics .................................................................................. 25
4.5
Setting mesh update parameters ..................................................................... 26
4.6
User defined function (UDF).......................................................................... 30
4.6.1
UDF for dynamic mesh .......................................................................... 30
4.7
Conclusion ...................................................................................................... 32
5
Modelling of rotary engine using dynamic mesh ................................................... 33
5.1
Dynamic mesh model with stationary rotor and deforming housing ............. 33
5.2
Dynamic mesh model with rotating rotor and stationary housing ................. 37
5.3
Eccentric rotation............................................................................................ 38
5.3.1
Rotation of a box .................................................................................... 38
5.3.2
Eccentric rotation of the box .................................................................. 40
5.3.3
Rotary engine with large housing........................................................... 43
5.4
Deforming housing......................................................................................... 46

iii

5.4.1
Butterfly Valve ....................................................................................... 46
5.4.2
Triangular valve...................................................................................... 50
5.4.3
Rotary engine with small clearance between rotor and housing ............ 54
5.4.4
Rotary engine with rotor and housing in perfect contact........................ 57
5.5
Conclusion ...................................................................................................... 61
6
Simulation of rotary engine .................................................................................... 62
6.1
Ideal gas compressible flow simulation.......................................................... 62
6.1.1
Rotary engine with small clearance........................................................ 63
6.2
Rotary engine with perfect contact................................................................. 66
6.3
Conclusion ...................................................................................................... 67
7
Conclusion .............................................................................................................. 68
8
Recommendation and Suggestion .......................................................................... 69
8.1
Recommendation for future work................................................................... 69
8.2
Suggestions for FLUENT............................................................................... 69
References ...................................................................................................................... 70
APPENDIX .................................................................................................................... 73
APPENDIX A - M-file used in the MATLAB for rotor and housing........................ 73
APPENDIX B - Triangular valve............................................................................... 78
APPENDIX C UDF for rotary engine..................................................................... 80

iv

TABLE OF FIGURES
Figure 2-1: Rotary Engine [1] .......................................................................................... 4
Figure 2-2: Peritrochoid curve of housing........................................................................ 5
Figure 2-3: Contour of the rotor ....................................................................................... 6
Figure 2-4: Parts of the rotor [1]....................................................................................... 7
Figure 2-5: Side housing [1]............................................................................................. 7
Figure 2-6: Out put shaft [1]............................................................................................. 8
Figure 2-7: Gas seal [1] .................................................................................................... 9
Figure 2-8: Intake port system [1] .................................................................................... 9
Figure 2-9: RENESIS (courtesy MAZDA) .................................................................... 10
Figure 2-10: Skycar (Courtesy Skycar) .......................................................................... 11
Figure 2-11: Micro-rotary engine used in Micro Electro Mechanical System (MEMS) 11
Figure 2-12: UAV (Courtesy UEL)................................................................................ 11
Figure 4-1: Spring based smoothing on interior nodes: Start [26] ................................. 20
Figure 4-2: Spring based smoothing on interior nodes: End [26] .................................. 20
Figure 4-3: Dynamic layering [26]................................................................................. 21
Figure 4-4: Results of splitting layer by constant height [26] ........................................ 22
Figure 4-5: Result of splitting layer by constant ratio [26] ............................................ 22
Figure 4-6: Remeshing at a deforming boundary [26] ................................................... 24
Figure 4-7: Dynamic mesh panel for smoothing ............................................................ 26
Figure 4-8: Effect of a spring constant factor of 1 on interior node motion [26]........... 27
Figure 4-9: Effect of a spring constant factor of 0 on interior node motion [26]........... 27
Figure 4-10: Dynamic mesh panel for layering.............................................................. 28
Figure 4-11: Dynamic mesh panel for local remeshing ................................................. 29
Figure 5-1: Complete RE (a) and a part of RE (b) that was considered......................... 33
Figure 5-2: Positions of the rotating rotor in a stationary housing [1]. .......................... 34
Figure 5-3: Positions of the housing with stationary rotor. ............................................ 37
Figure 5-4: Initial position of a square box. ................................................................... 38
Figure 5-5: Square box after 1 sec rotation. ................................................................... 39
Figure 5-6: Eccentric rotation of small box.................................................................... 41
Figure 5-7: Eccentric rotation of small box after 0.5 sec. .............................................. 41
Figure 5-8: Eccentric rotation of a small box after 1.0 sec............................................. 42
Figure 5-9: Eccentric rotation of a small box after 1.5 sec............................................. 42
Figure 5-10: Eccentric rotation of a small box after 2.0 sec........................................... 42
Figure 5-11: RE with large housing. .............................................................................. 43
Figure 5-12: RE with large housing after 0.5 sec. .......................................................... 44
Figure 5-13: RE with large housing after 1.0 sec. .......................................................... 44
Figure 5-14: RE with large housing after 1.5 sec. .......................................................... 44
Figure 5-15: RE with large housing after 2.0 sec. .......................................................... 45
Figure 5-16: RE with large housing after 2.5 sec. .......................................................... 45
Figure 5-17: RE with large housing after 3.0 sec. .......................................................... 45
Figure 5-18: Butterfly valve ........................................................................................... 48
Figure 5-19: Butterfly valve after 1.5 sec....................................................................... 49
Figure 5-20: Butterfly valve after 2.5 sec....................................................................... 49
Figure 5-21: Triangular valve......................................................................................... 50
Figure 5-22: Triangular valve after 2.5 sec with old macro. .......................................... 50

Figure 5-23: Triangular valve after 20 sec. .................................................................... 53

Figure 5-24: RE with small clearance ............................................................................ 54
Figure 5-25: RE with small clearance after 0.5 sec........................................................ 55
Figure 5-26: RE with small clearance after 1 sec........................................................... 55
Figure 5-27: RE with small clearance after 1.5 sec........................................................ 55
Figure 5-28: RE with small clearance after 3 sec........................................................... 56
Figure 5-29: Rotary engine having perfect contact between rotor and housing............. 57
Figure 5-30: Comparison of circle and rotor housing. ................................................... 58
Figure 5-31: RE with perfect contact at 0.2 sec. ............................................................ 58
Figure 5-32: RE with perfect contact at 0.5 sec. ............................................................ 59
Figure 5-33: RE with perfect contact at 1.0 sec. ............................................................ 59
Figure 5-34: RE with perfect contact at 1.5 sec. ............................................................ 59
Figure 5-35: RE with perfect contact at 2.4 sec. ............................................................ 60
Figure 6-1: Sketch of a rotary engine ............................................................................. 63
Figure 6-2: Static pressure simulation of rotary engine with clearance at 0.3 sec. ........ 64
Figure 6-3: Velocity vector simulation of rotary engine with clearance at 0.3 sec........ 64
Figure 6-4: Static pressure simulation of rotary engine with clearance at 0.5 sec. ........ 65
Figure 6-5: Velocity vector simulation of rotary engine with clearance at 0.5 sec........ 65
Figure 6-6: Static pressure simulation of rotary engine with perfect contact at 0.33 sec.
................................................................................................................................ 66
Figure 6-7: Velocity vector simulation of rotary engine with perfect contact at 0.33 sec.
................................................................................................................................ 66

vi

TABLE OF EQUATIONS
x = E cos(3 A) + ( R + C ) cos( A)
y = E sin(3 A) + ( R + C ) sin( A)

.Equation 2-1..................... 5

x = R cos(2V ) + (3E 2 / 2 R)(cos(8V ) cos(4V )) + E (1 9 E 2 / R 2 (sin 2 (3V )))1 / 2 (cos(5V ) + cos(V ))

y = R sin( 2V ) + (3E 2 / 2 R)(sin(8V ) sin( 4V )) + E (1 9 E 2 / R 2 (sin 2 (3V )))1 / 2 (sin(5V ) sin(V ))
Equation 2-2 .................... 6
x = R cos( 2V ) (3 E 2 / R )(sin( 6V ) sin( 2V )) + 2 E (1 ((9 E 2 / R ) sin 2 (3V ))) 1 / 2 cos( 3V ) cos( 2V )
y = R sin( 2V ) (3 E 2 / R )(sin( 6V ) cos( 2V )) + 2 E (1 ((9 E 2 / R ) sin 2 (3V ))) 1 / 2 cos( 3V ) cos( 2V )
. Equation 2-3 .................... 6
hmin > (1 + s )hideal
. Equation 4-1................... 21

hmin < c hideal

x = r cos( t )
y = r sin( t )

... Equation 4-2 .................. 22

Equation 5-1 .................. 40

x = r sin( t )

..Equation 5-2 .................. 40

y = r cos( t )
x = D cos

..Equation 5-3 .................. 46

y = D sin
x / D = cos

..Equation 5-4 .................. 47

= cos 1 ( x / D)
y = D sin(cos 1 ( x / D))

..Equation 5-5 .................. 47

vii

TABLE OF SYMBOLS AND NOTATIONS

C
A
V

omega
position
x, y, z
CG
CFD
time
E
d
R
hideal

c
s
vel
hmin
name
dt
r
D
rad/sec
RE
sec
3D
t
dtime
TDC
2D
UAV
UDF

Amount of parallel transfer

Angle of rotation for housing
Angle of rotation for rotor
Angular velocity
Angular velocity used in FLUENT
Array of coordinates
Cartesian coordinate of x, y and z axis respectively
Centre of gravity
Computational fluid dynamics
Current time supplied by FLUENT
Distance of eccentricity
Domain pointer
Generating radius
Ideal cell height
Layer collapse factor
Layer split factor
Linear velocity used in FLUENT
Minimum cell height
Name supplied by the user
Pointer to the dynamic zone data structure
Radius from the rotation
Radius of the circle
Radians per second
Rotary engine
Second
Three-dimensional
Time
Time step supplied by FLUENT
Top dead centre
Two-dimensional
Unmanned aerial vehicle
User defined function

viii

1 Introduction
Dr. Felix Wankel developed the first rotary engine in 1954 in NSU, Germany. Later it
has been licensed to many car makers like Ford, General Motors, Citroen, Mazda, etc.
Since last decade Mazda is the only car maker using rotary engine. Recently Mazdas
rotary engine, RENESIS was awarded International engine of the year 2003.
Rotary engine operation is similar to the typical reciprocating engine where four strokes
are completed by the rotation of the rotor inside a housing which is of complex
geometry called peritrochoid. Small size, low weight, high specific power density and
multi fuel capabilities are some of the advantages of rotary engine over reciprocating
engine. Further details of the rotary engine are discussed in chapter 2.
In order to design efficient rotary engine first it is necessary to understand how various
engine parameters affect the fuel-air mixing, fluid flow and combustion inside the
combustion chamber. Various mathematical models have been developed over the years
to enhance the understanding.
Automotive industry is using CFD as one of its research and design tools. Whether it is
the study of the external flow dynamics over the body of a vehicle or for the internal
flow through the engine/hood, CFD is helping automotive engineers to better
understand the physical flow processes and in turn to design improved vehicles.
FLUENT the most widely available commercial CFD software is used for the modelling
and simulation.
An extensive literature review has been carried out in chapter 3 to understand the
previous work done on rotary engine using CFD. Lack of experimental data which has
been acknowledged by all researchers hampers the proper validation of the model.
Some measured data which are available like flow pattern will be used for the
validation, once complete basic model is simulated.

It was found out that the last analysis of rotary engine using CFD was carried out in way
back in 1994. The shear development in CFD during this decade gives a lot of scope for
doing this project.
Dynamic capabilities of FLUENT and user defined function (UDF) are explained in
chapter 4. Various dynamic mesh update methods and their applicability were
discussed. User defined function is written in C programming language that can be
dynamically loaded with the FLUENT solver to enhance the standard features of the
code. UDF are defined using DEFINE macros that are supplied by FLUENT company.
Project Objective was

To setup a flexible model using dynamic mesh capabilities for rotary engine
First it was planned to set up a model having stationary rotor and rotating housing as it
was very simple to do. As this model was not flexible enough this model was not
pursued further.
Dynamic mesh option available in the FLUENT 6.1 was explored and UDF was used in
setting up the rotary engine model with stationary housing and rotating rotor.
The project was divided into many steps:
Eccentric rotation of the rotor: A rotary engine model is set up with oversized housing
having eccentric rotation of the rotor.
Deforming housing: Then the problem arising due to the contact between rotor and
housing is handled almost successfully.
Both the models with stationary rotor and rotating rotor are discussed in chapter 5.
Simulation is carried out with a model having no inlet and exhaust ports and rotor
pockets. The idea was to demonstrate simulation of rotary engine which is done in
chapter 6.
Based on the conclusion made in chapter 7, finally suggestion and recommendations are
made in chapter 8 for future work in setting up a practical rotary engine model.

2 Review of Rotary engine

Rotary engine is known for fewer mechanical parts and smooth operation. In one stroke
all the typical four stroke of an otto cycle is completed. Equations of rotor and rotor
housing and the basic constructions are explained in this chapter. This chapter ends with
the applications of rotary engine and scope for improvement.

2.1 History of Rotary engine

To overcome the shortcomings of around eighty years old reciprocating engines,
untiring research and development were carried out for rotary engine. Dr. Felix Wankel
along with the cooperation of NSU (West Germany) successfully developed a rotary
engine in 1954. Rotary engine is also called as Wankel engine. In spite of many
difficulties and problems it encountered first production of car having rotary engine was
in 1963, called NSU Spider. Later rotary engine was sub-licensed to General Motors,
Ford, Citroen, Mazda and other car makers. After many years of development Mazda
re-launched the rotary engine called RENESIS with the new RX-8 sports car. RENESIS
was awarded The international engine of the year 2003 and already meets EURO-4
emissions.

2.2 Principle of Operation

The operation of Rotary engine (RE) is completed by the typical four stroke similar to
the reciprocating engine but in one cycle. In the exhaust stroke, as the rotor rotates in
the clockwise direction the exhaust port opens and the residual gas moves out of the
combustion chamber due to the pressure difference. As the rotor rotates further in the
clockwise direction the intake port opens and fuel-air mixture moves into the
combustion. As the rotor turns further the intake port closes, the charge is compressed
and goes into the expansion stroke after being ignited near the top dead centre (TDC).
Most of the combustion is completed as the rotor moves past the TDC and later most of
the combustion products are driven out of the combustion chamber through the exhaust
port as the rotor rotates. The rotor turns eccentrically at one third of the crank shaft
speed.

Figure 2-1: Rotary Engine [1]

2.3 Salient features of rotary engine

The salient features of the RE over reciprocating engine are,
1. RE has no reciprocating parts: No problem of unbalance caused by the inertia of
reciprocating parts and engine vibration. Thus in RE vibration is very low as it is
possible to perfectly balance the engine.
2. RE requires no cranking mechanism: This in turn leads to smooth motion, less
mechanical loss, simple construction and compactness.
3. RE has no intake-exhaust valve mechanism: No mechanical noise generated by
the opening and closing of valves and also valves themselves are obstruction for
air flow. In RE opening and closing of ports are accomplished by the rotation of
the rotor.
Other features of RE are higher air flow capacity, light weight, compactness,
absence of connecting rod, flat torque characteristics, less moving parts and multifuel capability [2].
However, disadvantages of RE are
1. Longer combustion time
2. Higher heat losses

3. Incomplete flame propagation due to higher surface to volume ratio.

2.4 Basic Equation of the rotary engine

The equation of the rotor and housing is given by [1] and [10]
Housing
The basic curve forming the inner surface of the rotor housing is called peritrochoid.

The equation of peritrochoid for the RE in practical application can be expressed by the
coordinates of point P(x, y) as
x = E cos(3 A) + ( R + C ) cos( A)
y = E sin(3 A) + ( R + C ) sin( A)

Equation 2-1

Figure 2-2: Peritrochoid curve of housing

Where E is the distance of eccentricity

A is the angle of rotation
R is the generating radius and
C is the amount of parallel transfer
The equation given by both [1] and [10] are same for housing.
Rotor

The inner envelope of the peritrochoid is a basic curve forming the contour of the rotor.
The equation given in [1] is as follows

x = R cos(2V ) + (3E 2 / 2 R )(cos(8V ) cos(4V )) + E (1 9 E 2 / R 2 (sin 2 (3V )))1 / 2 (cos(5V ) + cos(V ))

y = R sin( 2V ) + (3E 2 / 2 R )(sin(8V ) sin( 4V )) + E (1 9 E 2 / R 2 (sin 2 (3V )))1 / 2 (sin(5V ) sin(V ))
Equation 2-2

The Equation 2-3 is given by [10]

x = R cos( 2V ) (3 E 2 / R )(sin( 6V ) sin( 2V )) + 2 E (1 ((9 E 2 / R ) sin 2 (3V ))) 1 / 2 cos( 3V ) cos( 2V )
y = R sin( 2V ) (3 E 2 / R )(sin( 6V ) cos( 2V )) + 2 E (1 ((9 E 2 / R ) sin 2 (3V ))) 1 / 2 cos( 3V ) cos( 2V )
. Equation 2-3

Figure 2-3: Contour of the rotor

Where E is the distance of eccentricity

R is the generating radius and
V is the angle of rotation which varies from
pi/6 to pi/2 for rotor face 1
5pi/6 to 7pi/6 for rotor face 2 and
3pi/2 to 11pi/6 for rotor face 3.
In figure 2-3 the rotor curve is generated using both equation 2-2 and 2-3. Equation 2-3
is preferred over equation 2-2 to generate rotor curve, as all the corners of the rotor
where in contact with the housing.

2.5 Basic construction of the rotary engine

The basic construction and parts unique to the RE are discussed in this section.

2.5.1 Rotor
The function of the rotor is similar to the
piston

and

reciprocating

connecting
engine.

rod

Rotor

of

directly

transmits the pressure of combustion gas to

the output shaft as rotating force. As the
rotor rotates, it automatically opens and
closes the intake-exhaust ports. A recess
provided on each rotor flank is called the
rotor recess. The compression ratio of an
engine is determined by its volume.
Figure 2-4: Parts of the rotor [1]

Rotor recess is also called as rotor pocket.

2.5.2 Housing
The housing of the RE is similar to the cylinder block and cylinder of the reciprocating
engine. As it is seen in figure 2-1, cooling air or
water is provided through the passages around
the outside of the inner surface of the rotor
housing. For higher rigidity and better cooling
effect, ribs or fins are also provided. Also in the
rotor housing spark plug holes, and exhaust port
(in the case of peripheral exhaust port) etc are
provided at their respective optimum positions.

Figure 2-5: Side housing [1]

The side housing is provided with an intake port (and exhaust port in RENESIS) as
shown in figure 2-5. The interior is made hallow for the passage of cooling water. For
rigidity and better cooling effect ribs are provided.

2.5.3 Phasing gear mechanism

For accurate control of the rotary motion of the rotor, RE is equipped with a phasing
gear mechanism. The peritrochoid profile motion of the rotor as discussed in the
section 2.4 is imparted by this mechanism.

2.5.4 Output shaft

The function of output shaft of a RE is similar to the crankshaft of the reciprocating
engine. It generates the rotating force by receiving the combustion pressure applied on
the rotor journal eccentric to the centre of rotation.

Figure 2-6: Out put shaft [1]

2.5.5 Gas sealing mechanism

The function of the gas sealing mechanism of a RE is similar to the piston and
compression rings of the reciprocating engine.
Gas seal mechanism consists of an apex seal to keep each working chamber gas-tight
from its adjacent ones, a side seal corresponding to the compression ring of the
reciprocating engine and a corner seal used at the junction of the above two seals. For
the peripheral port, it also plays the roles of both intake and exhaust valves.
Spring is provided on the back side of the each seal which will allow the seal to keep
close contact with the surface to be sealed even when it is worn.

Figure 2-7: Gas seal [1]

2.5.6 Intake and exhaust port system

Ports provided in the rotor housing are called peripheral port system. Ports provided in
the side housing are called side port system. The ports are opened and closed by the
apex seal.
In the peripheral intake system at high speed or heavy
loading due to the longer open time of the intake port
and the same direction of the intake flow and rotor
revolution contributes to a higher charging efficiency
and output. During low speed or light load combustion
rate is compromised due to the presence of greater
amount of burned gas caused by a longer overlap time
between the intake and exhaust ports.

Figure 2-8: Intake port system [1]

In the side intake system the direction of intake is normal to the rotor revolution thus
easily generates swirl and shorter overlap time between the intake and exhaust ports
results in stable combustion even during low-speed or light loading conditions. But at
high speed or heavy loading condition performance is lower than that of the peripheral
intake system.
Side exhaust port system will cause a great amount of leak of high temperature exhaust
gas into the rotor side space causing deterioration of the side and oil seals.

2.6 RENESIS

The key feature of the RENESIS engine

is side exhaust port system located at the
side housing as show in figure 2-9. The
chief advantage of this system is
elimination of the overlap between intake
and exhaust ports. Thus promotion more
stable

combustion

economy.

RENESIS

and
also

better

fuel

has

two

exhaust ports per rotor chamber thus

affording delay in the opening of the
exhaust ports which in turn leads to a
longer expansion cycle giving improved
Figure 2-9: RENESIS (courtesy MAZDA)

thermal efficiency, power output and fuel economy. RENESIS has two exhaust port per
rotor chamber thus able to afford delay in the opening of the exhaust port which in turn
leads to a longer expansion cycle giving improved thermal efficiency, power output and
fuel economy. Additionally, the intake port close timing has been ex tended, resulting in
increased charging volume and more power. RENESIS gains 30% in intake port area
over the previous engine, and this, combined with the delayed intake port close timing,
makes for a sizable increase in charging volume resulting in greater power output.
RENESIS employs a new cut-off seal located between the rotor's dual oil seals and side
seal. This sealing arrangement eliminates blow-by between intake and exhaust ports and
prevents carry-over of exhaust gas to the next intake cycle. Exhaust gas build-up against
the side seal can easily cause carbonization, but with the wedge-shaped or cuneiform
side seal, the seal shape is optimized to remove carbon.

2.7 Other applications of rotary engine

Many believe that RE has become extinct and few people know it is used in Mazda RX8. But the fact is application of RE can be found in many areas as discussed below.

10

Skycar

The Skycar volantor developed by Moller

International is capable of vertical take-off and
landing (VTOL) much as a helicopter and flies
from point of departure to destination much like
an airplane. Eight RE are used in the production
model. RE is used in Skycar due to high powerto-weight ratio at a reasonable cost and very
small for their power output.
Figure 2-10: Skycar (Courtesy Skycar)

Micro-rotary engine

A research has been carried out in University

of California at Berkeley, to develop a micro
RE fuelled by liquid hydrocarbons that will
last 7-14 times longer than the conventional
batteries (lithium or alkaline ). Also an engine
with 20% efficiency can be 10 times smaller
than the battery and still deliver the same
amount of energy.
Figure 2-11: Micro-rotary engine used in Micro Electro Mechanical System (MEMS)
(Courtesy University of California, Berkeley)

Unmanned Arial Vehicle (UAV)

UAV Engines Ltd (UEL) is using RE to power

Unmanned Arial Vehicle (UAV) because of the
compact size and quite running.
RE are also used as auxiliary power unit for
airplanes, high speed boat, industrial engines,
chain saws, etc.
Figure 2-12: UAV (Courtesy UEL)

11

2.8 Scope for improvement

The leakage caused by the wear of the seal at the corner of the triangular rotor has been
reduced by using ceramic seals. As mentioned earlier the emission problem has been
met recently using side exhaust port instead of peripheral port [3] in RENESIS.
RE typically consumes more fuel than the piston engine as the thermodynamic
efficiency of the engine is reduced by the long combustion chamber shape and low
compression ratio. Unique geometry of RE has the major advantage for operation using
multi fuels. In order to design efficient REs, it is necessary to have a good
understanding of how engine design and operating parameters affect the unsteady
multidimensional fluid flow, fuel-air mixing and combustion occurring inside the RE
combustion chamber.
Mathematical models have been developed in the last couple of decades to enhance
understanding of the flow field, fuel-air mixing and combustion in RE. Earlier
improvements in performance and emission have been achieved primarily through
experimental trial and error procedures making the analysis time consuming and costly.
In numerous design analyses in terms of rotor pocket geometry, spark plug position,
fuel injection can be carried out without much expense.

2.9 Conclusion
The principle operation and basic construction of RE has been discussed. Other than
automotive, RE are used in UAV, Skycar, etc. RENESIS really brought a dramatic
changes in the performance of RE. A literature review of the previous work done on RE
using CFD is done in the next chapter.

12

3 Review of the analysis of rotary engine using

computational fluid dynamics
Automotive industry is using Computational fluid dynamics as one of its research and
design tools. A literature review is carried out to take stock of the work done on RE
using CFD.

3.1 Computational Fluid Dynamics

Computational Fluid Dynamics (CFD) is the analysis of systems involving fluid flow,
heat transfer and associated phenomena such as chemical reaction by means of
computer based simulation.
Advantages of CFD over experimental analysis are,
1. Reduction of time and cost
2. Analysis of system under hazardous conditions which would have been difficult
experimentally
3. Unlimited results can be obtained.
CFD will never replace theoretical or experimental approach. It only helps to interpret
and understand the results of theory and experiment and vice versa [4].
Automotive industry is using CFD as one of its research and design tools. Whether it is
the study of the external flow dynamics over the body of a vehicle or for the internal
flow through the engine/hood, CFD is helping automotive engineers to better
understand the physical flow processes and in turn to design improved vehicles.
The determination of proper grid for the flow over or through a given geometric shape
is a serious matter. The way such a grid is determined is called grid generation. The
matter of grid generation is significant in CFD. Because the accuracy of the model is
entirely depends the type of grid one choose for a given problem [4].

3.2 Review of the analysis of RE using CFD

Predominantly the CFD aspect of the technology development program for predicting
the complex flow patterns occurring inside RE has been considered. Using CFD fuel

13

injection, spray and nozzle optimizations, rotor pocket and nozzle relocation, fluid flow,
fluid-air mixing and combustion can be analysed.
The first two dimensional model of the motored RE in the absence of combustion was
modelled by Shih et al. 1986 [5]. Their model is based on the computer code UF-LRC2D (University of Florida Lewis Research Center 2-Dimensional).
Grasso et al. [6] developed the first three dimensional model of the stratified charge RE
during early stages of flame propagation using REC-3D-FSC-86 computer code at John
Deere & Co., REC-3D-FSC-86, is a modified version of the KIVA code developed at
Los Alamos National Laboratory for the modelling of reciprocating engines [7].
Unfortunately detailed measurements are almost non-existent for REs [8]. In 1988 [9],
Abraham et al. compared the measured and computed chamber pressures thus for the
first time validating the three-dimensional computations for RE. They concluded that
the main feature of the combustion flow field in the engine is slow and non uniform
mixing of fuel and air that leads to long and incomplete combustion. And the cause for
non uniformity is attributed to the large fluid acceleration caused by the motion of the
rotor.
Using the computer code UF-LRC-3D which was an extension of UF-LRC-2D,
Steinthorsson et al [10] developed a model to study unsteady three-dimensional flow
field inside RE. They showed the velocity field inside a motored RE, the mixing of non
homogeneous fuel-air mixtures that enter trough the intake port, and the mixing that
takes place when the gaseous fuel is injected into the combustion chamber during
compression. But UF-LRC-3D code was unconditionally unstable for 3-dimensional
modelling.
Later Abraham and Bracco [11] carried out the comparisons of measured and computed
mean velocities and of measured fluctuating intensities and computed turbulence
intensities in a motored RE. The measurements were obtained by Laser Doppler
Velocimetry (LDV) at the Sandia National Laboratories [12, 13]. Results nearly agreed

14

with the experimental results and also they concluded that the turbulence intensity tends
to be rather homogeneous and of similar magnitude before and after TDC but
significantly inhomogeneous and of greater magnitude around TDC. This study also
supported their previous papers [7, 9].
A three-dimensional model for flow and combustion was applied to a direct-injection
stratified charge RE by Abraham and Bracco in 1989 [14] to identify the parameters that
control its burning rate. To enhance vaporization and the production of flammable
mixture, it is concluded that the orientation of the six sprays of the main injector with
respect to the air stream is important and no spray should be in the wake of any other
spray. It was predicted that if such a condition is respected, the indicated efficiency
would increase by some 6% at higher loads and 2% at lower loads. They further carried
out the work [19] on improvement in combustion with a pilot injector and predicted to
increase the indicated efficiency by about 5% over the range of operating conditions
studied.
Raju and Willis developed a new computer code for predicting the turbulent, and
chemically reacting flows with sprays occurring inside of a stratified-charge RE in 1990
[15]. One apparent finding of their study is that vaporization appears to be more rate-

controlling than mixing during the combustion process, in the case with advanced fuel
injection and spark ignition. Their results were in agreement with the experimental
flow-visualization pattern of Hamady et al. [16].
Comparison of computed and measured pressure in a premixed charge natural-gasfuelled RE, by Abraham and Bracco [23]. Computed and measured pressure agreed
within 12%. In addition they showed that the laminar conversion time becomes
dominant as flames approach walls, thus slowing down their propagation rate.
On their earlier developed computer model [6], Abraham and Bracco extended the
computations for premixed charge natural gas combustion of RE [17]. They concluded
that a leading pocket with a leading deep recess, a compression ratio of 9.5 and two
spark plugs, one before the minor axis and the other after, provide the best indicated

15

efficiency of the configurations. Also they mentioned that it is very likely that even
greater efficiency gains will be obtained if more drastic geometric changes are
investigated.
Abraham and Bracco [20], worked on ignition strategies in a stratified charge RE. They
predicted that a second ignition source located upstream of the main injector if ignites
the mixture consistently will lead to an increase in indicated efficiency of 6-8%. Also
they predicted that the gain in efficiency of the engine with two ignition sources would
be 7-10%, instead of 6-8% if a two-hole pilot injector is also used instead of the onehole pilot.
Again on the modified geometries of the rotor pocket, Abraham and Bracco along with
Epstein [21] studied two modifications to the pocket geometry of a stratified-charge RE.
In one of the modification, a pocket which is located towards the leading edge of the
rotor has been shown to produce recirculation within the pocket and in turn results in
faster burning of the fuel. In the second modification standard single pocket
configuration is replaced by two pockets with two injectors and two spark plugs. They
stated that such a configuration appears to result in better utilization of air in the
chamber and hence higher efficiency.
Finally in 1994 Abraham and Magi [18] carried out computations of the flows in the
cavity are presented for different configurations of the pilot injector and spark plug and
for different timings and fuel injection rates from pilot and main injectors. The results
showed that in order to maintain fast burning in the main chamber, the variables has to
be optimized so that the flame in the cavity ignites the charge in the main chamber prior
to TDC. Also, the relative timings of the main injector and pilot injector should be
optimized such that there is sufficient amount of flammable charge present in the
vicinity of the cavity opening by TDC.
Airflow visualization is carried by Hamady et al [16] in a motored stratified charge RE
for naturally aspirated and supercharged conditions. Major flow feature identified
include leakage past the apex seal, separation at the rotor pocket and a large

16

recirculation during intake. Also Hamady et al [22] carried out measurement in the
major plane near the central housing wall to quantify the blow-by (leakage at apex
seals) that has proven to be a dominating feature in the intake flow pattern. They
concluded that the intake stroke flow pattern was found to be dominated by a large scale
vertical structure circulating in the opposite direction from the rotor motion. The
compression stroke flow pattern was characterized by a rotor driven flow, although flow
reversal were present. These analyses are carried out to compare with a computer
model. They suggested simulation modelling of three combustion chambers, rather than
the single combustion chamber currently being used to include the effect of leakage at
apex seals. They observed leakage at apex seals even in modern RE.
Heywood et al [24] developed two types of computer models to analyze the
performance of both premixed-charge and direct-injection stratified-charge RE. In the
first model the combustion heat release rate as a function of crank angle is predicted
from measured engine chamber pressure data throughout the combustion process. The
second model predicts the development of the chamber pressure and hence engine
performance from an assumed mass burning rate or heat release profile.
Fuel-air mixing and combustion in a 2-D RE is carried out by Ramos et al [25]. The
effects of mixture stratification at the intake port and gaseous fuel injection on the flow
field and fuel-air mixing in a 2-D RE model have been investigated. They have shown
that the fuel distribution in the combustion chamber is a function of the air-fuel mixture
fluctuations at the intake port. The fuel is adverted by the flow field induced by the rotor
and is concentrated near the leading apex during the intake stroke. During compression,
the fuel concentration is highest near the trailing apex and lowest near the rotor. The
penetration of gaseous fuel injected into the combustion chamber during the
compression stroke increases with the injection velocity and results in recirculation
zones between the injector and trailing apex. And finally the fuel concentration near the
trailing apex and rotor is small except at high injection velocities.

3.3 Conclusion

17

The last computations of RE was carried out a decade from now (that is in 1994).
During this time a lot of development has happened in the CFD software field. This
itself provides an ample scope for the computation of RE.
As suggested by Hamday et al [22] analysis of RE must be carried with all the three
combustion chambers, rather than the single combustion chamber currently being used
to include the effect of leakage at apex seals.
KIVA code developed at Los Alamos National Laboratory or other codes used earlier
for the analysis of RE was not used because of the lack of supervision in using these
codes. Fluent 6.1, the most widely used CFD software was used for the modelling and
simulation of RE.

18

4 Dynamic mesh
FLUENT is generic CFD software. Dynamic mesh capability is the latest addition to
this software to make it more versatile. Applicability of the different dynamic mesh
updates methods is discussed. To enhance the standard features of the code user defined
functions are used.

4.1 Introduction to dynamic mesh

To model flows where the shape of the domain is changing with time due to motion in
the domain boundaries, dynamic mesh model can be used in FLUENT. The motion can
be a prescribed motion or unprescribed motion where the subsequent motion is
determined based on the solution at the current time. Based on the new position of the
boundaries at each time step the update of the volume mesh is handled automatically by
FLUENT. The motion can be described using either boundary profiles or user defined
functions (UDF) in FLUENT.

4.2 Dynamic mesh update methods

In FLUENT to update the volume mesh in the deforming region subjected to the motion
defined at the boundaries, three mesh motion methods are available.
1. Spring based smoothing
2. Dynamic layering and
3. Local remeshing

4.2.1 Spring based smoothing

The edges between two mesh nodes are idealized as a network of interconnected
springs. The equilibrium state of the mesh is referred to the initial spacing of the edges
before any boundary motion. Based on Hooks law at a given boundary node will
generate a force proportional to the displacement along all the springs connected to the
node.
The spring based smoothing for a cylindrical cell zone where one end of the cylinder is
moving is shown in figure 4-1 and 4-2.

19

Figure 4-1: Spring based smoothing on interior nodes: Start [26]

Figure 4-2: Spring based smoothing on interior nodes: End [26]

Applicability of the spring based smoothing method

Spring based smoothing can be used for any cell or face zone whose boundary is
moving or deforming.
For non-tetrahedral (non-triangular for 2D) when the following conditions are met,
spring based smoothing method is recommended:
1. The boundary of the cell zone predominantly in one direction. In other words
no anisotropic stretching or compression of the cell zone.
2. Predominantly the motion is normal to the boundary zone.

20

In any case the above mentioned conditions are not met, the resulting cell may have
high skewness value, since not all possible combinations of node pairs in nontetrahedral (non-triangular for 2D) are idealized as springs.

4.2.2 Dynamic layering

Dynamic layering can be used in prismatic (hexahedral and/or wedge) mesh zone to add
or remove layers of cells adjacent to a moving boundary, based on the height of the
layer adjacent to the moving surface. This method allows one to specify an ideal layer
height on each moving boundary. Based on the height h of the cells in layer j the
adjacent layer of cells is split or merged with the layer of cells next to it (layer i in
figure 4-3 ).

Figure 4-3: Dynamic layering [26]

The cell height are allowed to increase until

hmin > (1 + s )hideal

. Equation 4-1

Where hmin is the minimum cell height of cell layer j,

hideal is the ideal cell height and

s is the layer split factor

The cells are split when the above condition is met based on the specific layering
method
Constant height: The cells are split to create a layer of cells with constant height
hideal and a layer of cells of height h hideal .

21

Figure 4-4: Results of splitting layer by constant height [26]

Constant Ratio: In this method the cells are split such that locally, the ratio of
the new cell heights is exactly s everywhere.

Figure 4-5: Result of splitting layer by constant ratio [26]

In the case of compression, the cells in the layer j are being compressed until
hmin < c hideal

Equation 4-2

Where c is the layer collapse factor

The cells are merged when the above condition is met. In other words, the cells in layer
in j are merged with those in layer i.

22

Applicability of the dynamic layering method

Dynamic layering method can be used to merge or split cells adjacent to any moving
boundary provided the following conditions are met:
1. Even though the cell zone contain mixed cell shapes, all cells adjacent to the
moving face zone are either wedges or hexahedral (quadrilateral in 2D).
2. Except when sliding interfaces are used, the cell layer must be completely
bounded by one-sided face zones.
3. In the case of boundary face zone with two sided walls, the wall and wall
shadow pair must be split and use the coupled sliding pair interface option to
couple the adjacent cell zones.
4. In the case of model having periodic face zones in the cell zone where
dynamic layering is used only the serial version of the solver can be used.
However parallel solver for dynamic layering can be used, if the periodic
zones as periodic non-conformal interfaces.

4.2.3 Local remeshing method

Spring based smoothing method is normally used on zones with tetrahedral (triangular
for 2D) mesh. In the case of large boundary displacement compared to the local cell
sizes, the cell quality can deteriorate or the cell can become degenerate. This will lead to
convergence problem as the consequence of invalidated mesh like negative cell volume,
when the solution is updated to the next time step.
To overcome this problem, FLUENT locally remeshes cells which violate skewness or
the size criteria. The new cells thus formed satisfy the skewness and the size criteria,
mesh will be locally updated otherwise new cells are discarded.
Each cell is evaluated and marked for remeshing by FLUENT, if it meets one or more
of the following criteria:

23

1. It is small than the specified minimum size.

2. It is larger than the specified maximum size.
3. It has skewness greater than the specified maximum skewness.
Also, FLUENT allows to remesh triangular and linear faces on a deforming boundary.
Based on the moving and deforming loop of faces, FLUENT marks the deforming
boundary faces for remeshing.
Consider a tetrahedral mesh of a cylinder as shown in figure 4-6 whose bottom wall is
moving. On the deforming boundary, a single loop is generated at the bottom end of the
cylinder where the nodes are moving. Similar in approach as discussed in dynamic
layering method (section 4.2.2), the height of the faces connected to the nodes on the
loop is analyzed by FLUENT and later splits or merges the faces depending on the
specified ideal face height and split /merge factor.

Figure 4-6: Remeshing at a deforming boundary [26]

Applicability of the local remeshing method

Local remeshing method can be used only in zones that contain tetrahedral or triangular
cells.
While using local remeshing the faces on the deforming face zone can be remeshed only
if the following conditions are met:

24

1. The faces are triangular.

2. The faces to be remeshed are all adjacent to moving loop (In other words
moving nodes).
3. The faces are on the same face zone, and form a closed loop.
4. The faces are not part of a symmetry or conformal periodic boundary.

4.3 Volume mesh update procedure

Based on the methods described in section 4.2.1 4.2.3 the volume mesh is updated
automatically by FLUENT. Spring based smoothing and local remeshing methods will
be used for the boundaries of a tetrahedral (or triangular for 2D) cell zones that are
moving to update the volume mesh in this zone. For the zones containing hexahedral
and/or wedge cells dynamic layering method will be used.
By visiting adjacent cell zones FLUENT automatically determines which method to be
used and set appropriate flags for the volume mesh update method to be used. If the
motion is specified for the cell zone, FLUENT will visit all the neighbouring cell zones
and set the flag appropriately. FLUENT will visit only the adjacent cell zone if the
motion of a boundary zone is specified. No volume mesh update method will be applied
if a cell zone does not have any moving boundaries.

4.4 Solid-Body Kinematics

Solid-body kinematics is used if the motion is prescribed based on the position and
orientation of the centre of gravity of a moving object. For both cell and face zone this
is applicable.
The motion of the solid-body can be specified by the linear and angular velocity of the
centre of gravity. Velocities can be specified either as profiles or user defined function
(UDF).
In addition to the motion description the starting location of the centre of gravity and
orientation of the solid body must be specified. At every time step FLUENT
automatically updates the centre of gravity position and orientation.

25

4.5 Setting mesh update parameters

Smoothing

Enable smoothing option under mesh methods in the dynamic mesh panel to turn on the
spring based smoothing. The relevant parameters are as shown in figure 4-7.

Figure 4-7: Dynamic mesh panel for smoothing

Spring stiffness is controlled by adjusting the Spring Constant Factor between 0 and 1.
There is no damping on the spring if the value is 0 and the boundary node displacement
has more influence on the motion of the interior nodes. A default level of damping is
imposed on the interior node displacements if the value is 1.
The effect of spring constant is illustrated in figure 4-8 and 4-9. Degenerate cells are
created with spring constant factor of 1 as shown in figure 4-8. If the spring constant
factor is set to 0 the original mesh position is recovered as shown in figure 4-9.

26

Figure 4-8: Effect of a spring constant factor of 1 on interior node motion [26]

Figure 4-9: Effect of a spring constant factor of 0 on interior node motion [26]

If the model contains deforming boundary zones, Boundary Node Relaxation can be
used to control how the node positions on the deforming boundaries are updated. A
value of 0 prevents deforming boundary nodes from moving which is equal to turning
off smoothing on deforming boundary zones and a value of 1 indicates under-relaxation.
The solution can be controlled by using the values of Convergence Tolerance and
Number of Iteration.
Dynamic layering

Enable the layering option under mesh methods in the dynamic mesh panel to turn on
the dynamic layering.

27

A cell layer is split/merge is controlled by specifying Constant Ratio or Constant Height

under option. The Split Factor, s (Equation 4-1) and Collapse Factor, c (Equation 42) determines when a layer of cell that is next to a moving boundary is split or merged
with the adjacent cell layer, respectively.

Figure 4-10: Dynamic mesh panel for layering

Local remeshing

Enable the remeshing option under mesh methods in the dynamic mesh panel to turn on
the local remeshing.
Based on skewness, size and height of the adjacent moving zones, cells are remeshed
under local remeshing. Desired skewness of the mesh is specified in the Maximum Cell
Skewness. Cells having skewness above the maximum cell skewness value is marked
for remeshing. Based on the size criteria cells with volume above maximum cell volume
and below minimum cell volume are marked for remeshing. Based on skewness the
marking of cells is done at each time step whereas marking based on size criteria is
performed between the specified Size Remesh Interval since change in cell size is
typically small over one time step.

28

Figure 4-11: Dynamic mesh panel for local remeshing

FLUENT replaces the remeshed cells only if the quality of the remeshed cells has
improved. But this behaviour can be over ruled by disabling Must Improve Skewness
under option.
If the Sizing Function under option is enabled FLUENT marks cells based on the size
distribution generated by the sizing function in addition to the marking cells based on
the minimum and maximum volumes. The size distribution can be controlled by
specifying the Sizing Function Rate and Sizing Function Variation.
Sizing Function Variation is used to control how large or small an interior cell can be
with respect to its closest boundary cell. A value of zero indicates constant size
distribution away from the boundary. A value of 0.5 indicates an interior cell size can be
at most 1.5 times the size of the closest boundary cell. A value of -0.5 indicates an
interior cell size can be half the size of the closest boundary cell.
Sizing Function Rate is used to control how rapidly the cell size varies from the
boundary. This value can be specified between -1 and 1. A value of 0 indicates a linear
variation of cell size away from the boundary. If the value is positive, slower transition

29

from the boundary to the specified size function variation value. If the value is negative,
faster transition from the boundary to the specified size function variation value.
The resolution of the sizing function can be controlled by using Sizing Function
Resolution. The size of the background bins used to evaluate the size distribution with
respect to the minimum characteristic length of the current mesh is determined by the
sizing function resolution.

4.6 User defined function (UDF)

User defined function is written in C programming language that can be dynamically
loaded with the FLUENT solver to enhance the standard features of the code. UDF are
defined using DEFINE macros that are supplied by FLUENT company. Every UDF
contains the header file udf.h at the beginning of the source code file (#include udf.h),
which allows definitions for DEFINE macros and other fluent provided macros and
function to be included during the compilation process. Using UDF, FLUENT can be
customized to fit into particular modelling needs.

4.6.1 UDF for dynamic mesh

This document does not imply responsibility on the part of Fluent Inc. for the
accuracy or stability of solutions obtained using UDFs that are either user-generated
or provided by Fluent Inc. Support for current license holders will be limited to
guidance related to communication between a UDF and the FLUENT solver. The
above statements are taken from FLUENT 6.1 UDF manual. The explanation and
examples given for DEFINE macros given in the manual were limited. Bulk of the time
was consumed in understanding the function and working of macros. As it was stated
earlier dynamic mesh being still a new topic it did not helped the cause at all.
The motion of the dynamic zone has to be defined. If the zone is a rigid body a profile
or UDF can be used to define the motion of the rigid body. If the zone is a deforming
zone define the geometry and the parameters that control face remeshing if applicable.

30

If the zone is deforming and moving at the same time UDF can be used to define the
geometry and motion of the zone which change with time.
DEFINE macros that can be used to control the behaviour of the dynamic mesh are
1. DEFINE_CG_MOTION
2. DEIFNE_GEOM
3. DEFINE_GRID_MOTION

4.6.1.1 DEFINE_CG_MOTION

Using this macro, motion of a particular dynamic zone can be specified by providing
linear (vel) and angular velocity (omega) at every time step. FLUENT uses these
velocities to update the node positions on the dynamic zone.
Macro: DEFINE_CG_MOTION ( name, dt, vel, omega, time, dtime)
name is supplied by the user. This name is used while hooking the UDF file to
the dynamic zone.
dt is is a pointer to the dynamic zone data structure
vel is used to specify linear velocity by the user
omega is used to specify angular velocity by the user
time is the current time supplied by FLUENT
dtime is the time step supplied by FLUENT
dt, vel, omega, time, dtime are variables that are passed by the FLUENT solver to the
UDF. The linear and angular velocities are returned to FLUENT by overwriting the
arrays vel and omega respectively.

4.6.1.2 DEFINE_GEOM

This macro is used to define the geometry of the deforming zone. Through spring based
smoothing or after local remeshing the nodes are repositioned by calling
DEFINE_GEOM UDF when FLUENT updates a node on a deforming zone.
Macro: DEFINE_GEOM (name, d, dt, position)

31

name is supplied by the user. This name is used while hooking the UDF file to
the dynamic zone.
d is the domain pointer
dt is a pointer to the dynamic zone data structure
position is an array of coordinates
d, dt, position are variables that are passed by the FLUENT solver to the UDF. After
projection to the geometry defining zone the new position is returned to the FLUENT
by overwriting the position array.

4.6.1.3 DEFINE_GRID_MOTION

While updating the node position on a dynamic zone, FLUENT assumes that there is no
relative motion between the nodes. If there is a need to control the motion of each node
independently, this macro can be used. For example based on the deflection due to fluid
structure interaction the position of each node can be updated.
Macros: DEFINE_GRID_MOTION (name, d, dt, time, dtime)
The terms used in this macro are already explained in the previous two macros. d, dt,
time, dtime are variables that are passed by the FLUENT solver to the UDF.
DEFINE_GRID_MOTION is not used in setting up the model, but the remaining two
macros are discussed further as and when they are used.

4.7 Conclusion
All the methods in the dynamic mesh update and its applications are discussed. A lot of
effort was put to understand the working of UDF for dynamic mesh. UDF for dynamic
mesh are further explained in the next chapter as and when they are used.

32

5 Modelling of rotary engine using dynamic

mesh
The volume of the combustion chamber change with time due to the eccentric motion of
the rotor inside the stationary housing as discussed in section 2.2. To model flows
where the shape of the domain is changing with time due to motion in the domain
boundaries, dynamic mesh model is used as explained in chapter 4.

5.1 Dynamic mesh model with stationary rotor and deforming

housing
The project started with the simple idea of simulating the flow in RE, by keeping the
rotor stationary and housing rotating.
The advantages of this method are
1. The problem in setting up eccentric rotation of the rotor can be overcome as the
rotor is stationary.
2. Only one combustion chamber is simulated, hence computationally less
demanding
3. Simplicity in setting up grid as only one combustion chamber is considered as
shown in figure 5-1.

(a)

(b)

Figure 5-1: Complete RE (a) and a part of RE (b) that was considered.

33

4. In rotating rotor with a stationary housing model, at the point of contact

between rotor and housing FLUENT gives a lot of problem like negative cells
because of the deterioration of the mesh. This is one of the prominent reasons to
opt for the stationary rotor method.
Grid is generated for the initial rotor-housing position (as shown in figure 5-1(b)).
Using UDF the boundary of the housing is defined to deform with time. Thus obtaining
change in volume of the combustion chamber with time. UDF macro is written such that
the coordinates of the housing is made to change with respect to time.

Figure 5-2: Positions of the rotating rotor in a stationary housing [1].

In a practical RE the positions of the rotating rotor inside a stationary housing is shown
in figure 5-2. As explained in section 2.2 crankshaft completes three rotations when
rotor makes one complete rotation. Also it is clear that rotor rotates eccentrically with
the help of phasing gear mechanism so that all the apex seal is in touch with the
housing. In doing so rotor is moving upward in figure 5-2(ii) as the rotor rotates in
clockwise direction in figure 5-2(iii) rotor is moving towards the right hand side.
Further in figure 5-2(iv) rotor is moving down wards and finally to the left hand side
(figure 5-2(vi)) as rotor rotates. This sequence is periodic with the rotation of the rotor.

34

To generate coordinates with stationary rotor and deforming housing model, housing is
made to rotate in the same sequence but in opposite direction. MATLAB is used to
rotate the housing as said in the above sequence using m-files. M-files used to generate
curves for rotor and housing is provided in APPENDIX A. Positions of the housing at
an interval of 30 degree rotation with stationary rotor are shown in figure 5-3. In the
figure 5-3 only the part of the housing coming in front of the considered phase of the
rotor is taken into consideration. These coordinate later could have been used in UDF to
define the boundary of the housing changing with time. Thus the volume of the
combustion chamber is made to vary with time.

(a): Zero degree rotation (Initial position)

(b): 30 degree rotation

(c): 60 degree rotation

(d): 90 degree rotation

35

(e): 120 degree rotation

(f): 150 degree rotation

(g): 180 degree rotation

(h): 210 degree rotation

(i): 240 degree rotation

(j): 270 degree rotation

36

(k): 300 degree rotation

(l): 330 degree rotation

Figure 5-3: Positions of the housing with stationary rotor.

This method was not pursued further due to the following limitations:
1. At the later stage inlet and exhaust ports can not be provided as this will make
the housing boundary discontinuous.
2. Ignition, inlet and exhaust ports effect can not be considered.
3. Also, the effect on the flow when the rotor rotates in the stationary housing is
not same as the effect of stationary rotor and rotating housing. This effect was
over looked because of the simplicity it offered.
According to Sierens et al. [27], using an electric discharge system while rotating the
housing around a fixed rotor, Wolters (1978) measured gas velocities in the combustion
chamber of a RE. No additional information about Wolters measurements is provided.

5.2 Dynamic mesh model with rotating rotor and stationary

housing
It was necessary to develop a basic model which handles almost all the complexity of a
practical RE. Again dynamic mesh model is used as the shape of the domain is
changing with time due to motion in the domain boundaries.
The setting up of the basic model is divided into three parts:
1. To set up a model having eccentric rotation of rotor.
2. To set up a model having a rotating part being in perfect contact with the
housing.

37

3. After successfully developing the above two model the same idea is
implemented into the RE model and carryout flow simulation.
Once the basic model is set up then it can be further developed like providing inlet and
exhaust port for the flow analysis, providing sparkplug for combustion analysis, varying
the rotor recess geometry, spark plug location, etc.

5.3 Eccentric rotation

First stage is to set up a model with only eccentric rotation of the rotor inside an
oversized housing. Again this stage is further divided
1. A model having a part rotating over its centre of gravity.
2. A model having a part rotating at a distance from the centre and also rotating
over its centre of gravity.
3. After the successful run of the above model, the same UDF is customized for the
RE model.

5.3.1 Rotation of a box

A square box is constructed at the centre of a big box as shown in figure 5-4. Rotation
of square box over its centre of gravity (CG) is defined using UDF.

Figure 5-4: Initial position of a square box.

38

Grid is generated using GAMBIT which is pre processor for generation of grid and
applying boundary conditions.
Using DEFINE_CG_MOTION macro, motion of a particular dynamic zone can be
specified by providing with the linear (vel) and angular velocity (omega) at every time
step. FLUENT uses these velocities to update the node positions on the dynamic zone.
#include "udf.h"
DEFINE_CG_MOTION(rotate, dt, vel,omega, time, dtime)
{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
}

The standard DEFINE_CG_MOTION features are discussed in section 4.6.1.1. NV_S

is a vector macros used to initialize the linear (vel) and angular (omega) velocities.
In this case CG motion of one rad/sec is specified using only angular velocity.
omega [2] = 1.0
The axis of rotation is specified in the square bracket. 0, 1 and 2 represents x, y and z
axis respectively.
Local remeshing and smoothing is used as the mesh is triangular but not much
concentration was given on the parameters of remeshing and smoothing. Figure 5-5
shows the preview of the square box after 1 sec. The quality of remeshing is not so good

Figure 5-5: Square box after 1 sec rotation.

39

5.3.2 Eccentric rotation of the box

Similar to the eccentric rotation of the rotor, in this model the box is made to rotate at a
distance from the centre of a big box by specifying angular velocity (omega). Also the
small box is made to rotate about its CG at three times the velocity of the eccentric
rotation as the rotor rotates three times during one complete cycle in RE. Using linear
velocity equation, rotation of the small box over its CG is specified.
A body rotating with angular velocity at a radius r from the rotation centre at
position r(t) has Cartesian coordinates
x = r cos(t )
y = r sin( t )

Equation 5-1

To find the velocity above equation is differentiated with respect to time

.

x = r sin( t )

..Equation 5-2

y = r cos( t )
In this case angular velocity is given three times the value of omega.
#include "udf.h"
#define E 3.0 /*Eccentricity*/
real beta;
DEFINE_CG_MOTION(eccentric, dt, vel, omega, time, dtime)
{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
beta = 3;
vel[0] = -beta*E*sin(beta*time);
vel[1] = beta*E*cos(beta*time);
vel[2] = 0;
}

In this case linear velocity (vel) is used to specify motion of the small box at a distance
E from the centre of big box.
vel[0], vel[1] and vel[2] represents x, y and z Cartesian coordinates respectively.
beta is the angular velocity which is assigned as three rad/sec.
As explained earlier angular velocity (omega) is used to specify the motion of the box
about its own CG.
The initial grid is generated using GAMBIT with a small box placed at a distance from
the centre of the big box as shown in the figure 5-6. Local remeshing and smoothing is

40

used as the mesh is triangular and again not much concentration was given on the
parameters of smoothing and remeshing. The preview position of the small box is
shown from figure 5-16 to 5-20m, fine mesh is used and quality of remeshing is not
good (white patches) as the small box rotates.

Figure 5-6: Eccentric rotation of small box.

Figure 5-7: Eccentric rotation of small box after 0.5 sec.

41

Figure 5-8: Eccentric rotation of a small box after 1.0 sec.

Figure 5-9: Eccentric rotation of a small box after 1.5 sec.

Figure 5-10: Eccentric rotation of a small box after 2.0 sec.

42

5.3.3 Rotary engine with large housing

After the successful modelling of the eccentric rotation of a small box, the same UDF is
customized for RE.
#include "udf.h"
#define E -0.01542 /* eccentricity of the rotor*/
real beta;
DEFINE_CG_MOTION(eccentric, dt, vel, omega, time, dtime)
{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
beta = 3;
vel[0] = -beta*E*sin(beta*time);
vel[1] = beta*E*cos(beta*time);
vel[2] = 0;
}

In the above UDF only the value of eccentricity E is changed.

In this case RE with large clearance between rotor and housing is considered to avoid
the problem of rotor and housing contact as shown in figure 5-11. Again not much
concentration was paid on the local remeshing and smoothing parameters. Finer mesh is
used and the remesh quality is not so good.

Figure 5-11: RE with large housing.

43

Figure 5-12: RE with large housing after 0.5 sec.

Figure 5-13: RE with large housing after 1.0 sec.

Figure 5-14: RE with large housing after 1.5 sec.

44

Figure 5-15: RE with large housing after 2.0 sec.

Figure 5-16: RE with large housing after 2.5 sec.

Figure 5-17: RE with large housing after 3.0 sec.

45

5.4 Deforming housing

After the successful modelling of eccentric rotation of the rotor, the problem arising
from the contact between the rotor and housing is taken in this section. As it was stated
earlier rotating rotor with a stationary housing model, at the point of contact between
rotor and housing FLUENT gives a lot of problem like negative cells because of the
deterioration of the mesh. If the clearance is small or in contact between the rotating
part and stationary housing the boundary nodes of the housing are pulled by the rotating
part. Thus the housing loses its original geometry and negatives cells are formed.
To overcome this problem a new macro, DEFINE_GEOM macro is used to define the
geometry of the deforming zone. FLUENT updates the node on the deforming zone
through spring smoothing or after local remeshing. The node is repositioned by calling
DEFINE_GEOM macro at each time step individually.
The problem is divided into different stages
1. First with a model having small clearance between rotating part and housing.
2. A model having rotating part and housing in perfect contact.
3. RE with small clearance between rotor and housing.
4. RE with rotor and housing being in perfect contact.

5.4.1 Butterfly Valve

A butterfly valve is considered with a small clearance between rotor and housing as
shown in figure 5-18. DEFINE_GEOM is used to define the geometry of the deforming
housing.
As the housing is perfect circle the equation of the circle is used in DEFINE_GEOM to
define the nodes of the deforming housing.
For a circle
x = D cos
y = D sin

..Equation 5-3

Where D is the radius of the circle

46

 is the angle of rotation

This is substituted as
x / D = cos

..Equation 5-4

= cos 1 ( x / D)
y = D sin(cos 1 ( x / D))

..Equation 5-5

The above equation is used to define the circular geometry of the housing.
#include "udf.h"
#define D 10 /* radius of the circle */
DEFINE_GEOM(slippery_arcs, domain, dt, position)
{
real y;
y=position[1];
if (y>=0.0)
position[1] = D*sin(acos(position[0]/D));
else
position[1] = -D*sin(acos(position[0]/D));
}
DEFINE_CG_MOTION(rotate, dt, vel,omega, time, dtime)
{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
}

DEFINE_GEOM is used along with DEFINE_CG to define the rotation of the valve.
The standard features of UDF DEFINE_GEOM are discussed in section 4.6.1.2.
Here position [0] and position [1] corresponds to x and y coordinates of the housing.
Due to the trigonometric limitation y coordinates less than zero are just multiplied by
negative sign.
Specified parameters are
Spring based smoothing
Spring constant factor

0.02

Boundary node relaxation

0.30

Convergence tolerance

0.001

Number of iteration

Local remeshing

47

Minimum cell volume (m3)

0.017

Maximum cell volume (m ) 0.230

Maximum cell skewness

0.65

Size remesh interval

Size function resolution

Size function variation

0.5

Size function rate

-0.5

Dynamic zone

Valve
Type

rigid body

Cell height

0.005

UDF file

rotate

Centre of gravity location

x=-0.0; y=0.0

Housing
Type

deforming

Geometry UDF

Slippery_arcs

Cell height

0.005

Height factor

0.5

Maximum skewness

0.65

The preview of the butterfly valve is show from figure 5-18 to 5-20.

Figure 5-18: Butterfly valve

48

Figure 5-19: Butterfly valve after 1.5 sec.

Figure 5-20: Butterfly valve after 2.5 sec.

At the end of 2.5 sec the mesh quality was reasonable good and also it was able to retain
the geometry of the housing in good shape.

49

5.4.2 Triangular valve

Triangular valve model is somewhat similar to the condition of RE (At any time rotor is
in contact with the housing at the three apex sealing with the housing) is constructed as
shown in figure 5-22.
As the housing is again a perfect circle the same macros used in butterfly valve is used
again.

Figure 5-21: Triangular valve.

Figure 5-22: Triangular valve after 2.5 sec with old macro.

50

Even after playing with local remeshing and spring based smoothing parameters, the
geometry of the housing was not able to retain to its original shape as shown in figure 522.
If we revisit the equations used DEFINE_GEOM, using the x-coordinates, ycoordinates of the nodes are repositioned. But as the node is displaced during the time
step both x and y coordinates may have moved away from the original circle
coordinates. So the x-coordinates used to find y-coordinates may itself be displaced
from the original circle coordinates thus giving a wrong y-coordinate value.
Also for housing having complex equation as in the case of RE, it is difficult to inverse
the equation (like Equation 5-4) and use this method.
To overcome this problem, x and y-coordinates for the nodes placed on the housing at
initial conditions are generated. At each time step, as the FLUENT calls
DEFINE_GEOM to reposition each node individually, the distance between the
displaced node and the generated points is compared to identify the closest and thus
reposition the node with the coordinates of the closest point. By doing so both x and y
coordinates are replaced instead of doing only one coordinates.
DEFINE_GEOM is used along with DEFINE_CG to define the rotation of the valve.
The standard features of UDF DEFINE_GEOM are discussed in section 4.6.1.2.
#include "udf.h"
#define NCOORDS 126 /* number of coordinates points*/
DEFINE_GEOM(slippery_arcs, domain, dt, position)
{
int i;
real dist2;
int closest_i;
real distmin =100.00;
real nodex[NCOORDS] = { ..... x coordinates of the housing
};
real nodey[NCOORDS] = {
};

.....y coordinates of the housing

for(i=0; i < NCOORDS; i++)

{
dist2 = (position[0] - nodex[i])*(position[0] - nodex[i])
+ (position[1] - nodey[i])*(position[1] - nodey[i]);
if (dist2 < distmin)
{
closest_i = i;

51

distmin = dist2;

}
position[0] = nodex[closest_i];
position[1] = nodey[closest_i];

DEFINE_CG_MOTION(rotate, dt, vel,omega, time, dtime)

{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
}

In the above programme nodex and nodey are generated x and y coordinates which are
defined as arrays.
NCOORDS is the number of coordinates points generated
i is the integer used for counting number of arrays.
dist2 is the square of the distance between the generated points and the displaced node
point.
distmin is the variable to store the minimum value of the square of the distance
between the generated points and the displaced node point.
When distmin reaches minimum value the corresponding x and y coordinates are used
to reposition the nodes on the housing using position.
Again position [0] and position [1] corresponds to x and y coordinates of each node on
the housing.
This same UDF along with x and y coordinates are given in APPENDIX B. Triangular
Valve after 20 sec is shown in figure 5-24.

52

Figure 5-23: Triangular valve after 20 sec.

Parameters

Spring based smoothing

Spring constant factor

0.02

Boundary node relaxation

0.30

Convergence tolerance

0.001

Number of iteration

Local remeshing
Minimum cell volume (m3)

0.017

Maximum cell volume (m3) 0.230

Maximum cell skewness

0.65

Size remesh interval

Dynamic zone

Valve
Type

rigid body

Cell height

0.4

UDF file

rotate

Centre of gravity location

x=-0.0; y=0.0

Housing
Type

deforming

Geometry UDF

Slippery_arcs

Cell height

0.4

53

Height factor

0.5

Maximum skewness

0.65

The quality of the remesh at the end of 20 sec was good and the geometry of the
housing was also in perfect shape.

5.4.3 Rotary engine with small clearance between rotor and housing

Figure 5-24: RE with small clearance

After the successful run of the triangular valve, RE with small clearance between rotor
and housing is modelled using GAMBIT as shown in figure 5-24. The UDF used in
triangular valve is implemented for both RE with clearance and perfect contact. Exactly
same UDF is used for both RE model. As the geometry of the housing is not as simple
as a circle, UDF is bit modified for the reason discussed in section in 5.4.4. UDF is
given in APPENDIX C.

54

Figure 5-25: RE with small clearance after 0.5 sec.

Figure 5-26: RE with small clearance after 1 sec.

Figure 5-27: RE with small clearance after 1.5 sec.

55

Figure 5-28: RE with small clearance after 3 sec.

Parameters

Spring based smoothing

Spring constant factor

0.5

Boundary node relaxation

0.4

Convergence tolerance

0.001

Number of iteration

Local remeshing
Minimum cell volume (m3)

3.088e-06

Maximum cell volume (m3) 1.577e-05

Maximum cell skewness

0.65

Size remesh interval

Dynamic zone

Rotor
Type

rigid body

Cell height

0.0039

UDF file

eccentric

Centre of gravity location

x=-0.01542; y=0.0

Housing
Type

deforming

Geometry UDF

Slippery_arcs

56

Cell height

0.0039

Height factor

0.5

Maximum skewness

0.65

The model is running properly even after more than 3 sec. Simulation is carried out on
this model first before doing on the RE with rotor and housing in perfect contact. This
gives an idea in setting up the parameters for simulation in RE with perfect contact.
Also the simulation in model having perfect will be slower than the RE with small
clearance.

5.4.4 Rotary engine with rotor and housing in perfect contact

Figure 5-29: Rotary engine having perfect contact between rotor and housing.

After the successful run of the triangular valve model, RE with perfect contact between
rotor and housing is modelled as shown in figure 5-29. The same UDF used for the
triangular valve is used with some modification for both RE with and without clearance
between rotor and housing model.
During grid generation nodes are placed at equal distance on the circumference of the
housing and using MATLAB coordinates are generated for equal increments of angle.
Increment of angle is set equal to the circumferential distance between the nodes. In the
case of triangular valve the housing was a perfect circle so nodes are uniformly placed
over the housing in terms of angle between the nodes from the origin.

57

Figure 5-30: Comparison of circle and rotor housing.

Where as in the case of rotor housing, the radius from the centre of the origin is not
constant as shown in the figure 5-31, thus the nodes placed at equal distance on the rotor
housing will not be at equal interval of angle from the origin. This in turn leads to the
mismatch in the coordinates generates at equal increment of angle using MATLAB and
the coordinates of the nodes placed during grid generation. To overcome this problem
coordinates are generated at closer interval than the actual node position. In other words
if the housing has 100 nodes, say 200 coordinates at equal interval of angle from the
origin is generated and fed into the UDF. The UDF is given in APPENDIX C.

Figure 5-31: RE with perfect contact at 0.2 sec.

58

Figure 5-32: RE with perfect contact at 0.5 sec.

Figure 5-33: RE with perfect contact at 1.0 sec.

Figure 5-34: RE with perfect contact at 1.5 sec.

59

Figure 5-35: RE with perfect contact at 2.4 sec.

Parameters

Spring based smoothing

Spring constant factor

0.2

Boundary node relaxation

0.5

Convergence tolerance

0.001

Number of iteration

Local remeshing
Minimum cell volume (m3)

3.053e-06

Maximum cell volume (m3) 1.483e-05

Maximum cell skewness

0.65

Size remesh interval

Dynamic zone

Rotor
Type

rigid body

Cell height

0.0039

UDF file

eccentric

Centre of gravity location

x=-0.01542; y=0.0

Housing
Type

deforming

Geometry UDF

Slippery_arcs

Cell height

0.0039

60

Height factor

0.5

Maximum skewness

0.65

Figure 5-31 to 35 shows the RE with zero clearance at different time position. At the
end of the preview the shape of the housing is retained. As it is said earlier the
coordinates generated is more the no of actual nodes on the housing, the quality of the
mesh near the housing was not as good as it was in the triangular valve model. In turn
preview failed at 2.4 sec due to negative cells. The quality of remeshing other than near
rotor housing is good.
Before opting for this model many other options like with finer mesh, coarser mesh,
varying intensity of the coordinates only at some region, increasing area of contact
between rotor and housing, etc has been checked out without any success. Also, the
coordinates of the node on the housing was tried to measure indirectly by placing vertex
on each node (In GAMBIT there is no option to check the coordinate value of the node)
and thus generating coordinate value. Due to compromise done while placing the vertex
finally the model with UDF having measured coordinates did not worked out at all.
After exploring many options it is settled for present UDF having generated coordinates
points double the number of nodes on the housing. The logic has already been tested
convincingly in triangular valve model. Only a fine tuning is required in terms of
coordinates generated which is used in UDF for smooth running.

5.5 Conclusion
Dynamic capabilities has been explored deeply to setup a flexible RE model which
accept complexities of practical RE like inlet and exhaust ports, ignition, etc., The
eccentric rotation of the rotor is accomplished. The problem arising from the contact
between the rotating part and the housing has been handled successfully as shown in
triangular valve. Finally a RE model with perfect contact between rotor and housing is
set up. The UDF used in triangular valve model was customized for RE model due to
the complex shape of the rotor housing. For demonstration simulation of the available
RE model has been carried which will be discussed in the next chapter.

61

6 Simulation of rotary engine

A one-dimensional model of combustion in premixed-charge RE was introduced in the
early seventies just after a corresponding model for reciprocating engine. But the former
was then extended to sprays and spray combustion well ahead of similar upgrading for
reciprocating engines. Since then, of course, modelling for reciprocating engines has
blossomed and for RE has withered [6].
The last published work on computation of RE was in 1994. Since last decade a lot of
development has happened in CFD. Initial objective will be to set up a flexible 2-D
model and validate the results with the experimental flow-visualization pattern of
Hamady et al. [16].
Setting up a model with an eccentric rotation of the rotor and handling the deforming
housing due to the contact between rotor and housing was a difficult phase as discussed
in chapter 5. Simulation was carried out with the model having no inlet and exhaust
ports and rotor pocket for demonstration. Hence only the effect of the rotor motion on
the flow pattern is observed.

6.1 Ideal gas compressible flow simulation

Simulation for Ideal gas compressible flow is carried out for both RE with small
clearance and perfect contact. Standard k turbulence model is used for its
robustness and accuracy for a wide range of turbulent flows. The simulations are carried
out for demonstration only as there were no inlet and exhaust ports and rotor pockets in
the model analysis can not be carried out with respect to practical RE. Hence much
importance was not given while assigning the parameters for simulation.
The parameters set for both RE with clearance and perfect contact are as follows:
Operating pressure 101325 pascal
Initialization
Gauge pressure

X Velocity

0.2 m/s

62

Y velocity

0.2m/s

Temperature

300k

Figure 6-1: Sketch of a rotary engine

During the simulation rotor is rotating in anti-clockwise direction as shown in figure 61.

6.1.1 Rotary engine with small clearance

In RE with small clearance the cells between the rotor and housing will not be
comparatively strained much when the rotor and housing comes closer. The iteration
parameters are as follows
Time step

1e-6

Maximum number of iteration per time step 20

The simulation of RE with small clearance are shown from figure 6-2 to 6-5. Leakage at
the apex sealing can be observed at all the time. Also velocity of the flow is maximum
at the sealing. Pressure started to build up due to compression in figure 6-2, but in figure
6-4 pressure reduces to minimum due to leakages at apex sealing.

63

Figure 6-2: Static pressure simulation of rotary engine with clearance at 0.3 sec.

Figure 6-3: Velocity vector simulation of rotary engine with clearance at 0.3 sec.

64

Figure 6-4: Static pressure simulation of rotary engine with clearance at 0.5 sec.

Figure 6-5: Velocity vector simulation of rotary engine with clearance at 0.5 sec.

65

6.2 Rotary engine with perfect contact

Figure 6-6: Static pressure simulation of rotary engine with perfect contact at 0.33 sec.

Figure 6-7: Velocity vector simulation of rotary engine with perfect contact at 0.33 sec.

66

After getting an idea about the time step input simulation of the RE model with perfect
time contact between rotor and housing is carried out. The iteration parameters are as
follows
Time step

1e-6

Maximum number of iteration per time step 100

As expected, iteration was slow in this case. Figure 6-6 and 6-7 shows the simulation of
rotary engine with perfect contact.
As shown in figure 6-6 pressure build up can be seen due to compression and no
leakage can be observed in figure 6-7.
Nothing serious analysis can be carried out with the available model as there were no
inlet and exhaust port and rotor pocket. But definitely this is first step in setting up a 2D
model with rotor pocket and inlet and exhaust port considering all the three combustion
chamber and validating the flow pattern with that of the experimental results [16].

6.3 Conclusion
Simulation of RE with small clearance and perfect contact between rotor and housing is
carried out for demonstration. Nothing really can be drawn in terms of flow analysis as
there were no inlet and exhaust port and rotor pocket. This is the first step in setting up a
2D model of a practical RE.

67

7 Conclusion
Principle operation of RE is completed by the typical four stroke of an otto cycle.
During this time the crankshaft rotates at three times the speed of the rotor. RE is
mechanically simpler IC engine with no valves, cranking mechanism and reciprocating
part. This makes the RE compact, smooth and powerful.
Review of the analysis of the RE using CFD is carried out, and lack of experimental
data is acknowledged. The last computational analysis on RE was carried out in 1994.
In order to design efficient RE it is necessary to under the influence of the engine
parameters on fluid flow, fuel-air mixing and combustion of RE. Dynamic mesh
capabilities and UDF of FLUENT are explored in setting up a flexible RE model.
Modelling of RE was started with stationary rotor and rotating housing. Later
modelling of RE with stationary housing and rotating rotor is set up as it was found
more flexible in handling inlet and exhaust ports and ignition for further work. The
eccentric rotation of the rotor is modelled successfully. The most challenging part of
this project was the problem arising from the contact between the rotating part and the
housing as in triangular valve is modelled to run smoothly as a result of tenacious effort.
Further the complex geometry of the rotor housing does not make the situation any
simple. Thus modelling of the RE was done with limited success.
Ideal gas compressible flow simulation is carried out both on RE model having small
clearance and perfect contact between rotor and housing. This was carried out to
demonstration the simulation of RE.

68

8 Recommendation and Suggestion

Based on the experience and knowledge gained during this project following
recommendations and suggestions are made to obtain better results.

8.1 Recommendation for future work

1. Perseverance will prove a key in modelling a smooth running RE using dynamic
mesh capabilities of FLUENT. The area needs to be looked after is the
generation of coordinates used in UDF.
2. With standard rotor pocket configuration, taking inlet and exhaust port into
consideration, 2D flow analysis has to be carried out first without combustion
3. Simulation modelling of three combustion chambers, rather than single
combustion chamber currently being used to include the effect of leakage at
apex seals as suggested by Hamady and validated with experimental results.
4. For a better analysis of RE 3D analysis must be carried out.

8.2 Suggestions for FLUENT

1. Stability of the dynamic mesh capabilities has to be improved.

During segmentation violation in mesh preview, FLUENT has to be restarted

again.

Sometimes when files are read from recent files, UDF goes missing in
dynamic zone.

Sometimes at the initial stage of the mesh preview, graphics window

becomes blank.

Couple of times case file having same parameters has resulted in different
mesh preview.

2. Understanding the UDF for dynamic mesh capabilities proved difficult.

3. Coordinates of the nodes can not be measured in the GAMBIT. Again GAMBIT
is not stable enough when compared to other CAD software.

69

References
1. Yamamoto, K., Rotary Engine, Sankaido, Tokyo, 1981.
2. Jones, C., Ellis, D.R and Meng, P. R., Multi-Fuel Rotary Engine for General
Aviation Aircraft, AIAA/SAE/ASME 19th Joint Propulsion Conference, Seatle, Wash.,
Paper No. AIAA-83-1340, 1983.
3. Ritsuharu Shimizu, Harno Okimoto, Seijo Tashima and Suguru Fuse, The
Characteristic of Fuel Consumption and Exhaust Emission of the Side Exhaust port
Rotary Engine, SAE International Congress and Exposition, Detroit, Mich., SAE Paper
No. 950454, 1995.
4. John D. Anderson, Computational fluid dynamics: The basics with applications,
McGraw-Hill, Inc. New York, 1995.
5. Shih, T.I-P., Yang, S-L., and Schock, H.J., A Two-dimensional Numerical Study of
the Flow Inside the Combustion Chamber of a Motored Rotary Engine, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 860615, 1986
6. Grasso, F., Wey, M.J., Bracco, F.V., and Abraham, J., Three-Dimensional
Computations of Flows in a Stratified-Charge Rotary Engine, SAE International
Congress and Exposition, Detroit, Mich., SAE Paper No. 870409, 1987.
7. Amsden, A.A., Ramshaw, J.D., ORourke, P.J., and Dukowicz, J.K., KIVA: A
Computer Program for Two and Three Dimensional Fluid Flows with Chemical
Reactions and Fuel Sprays, LA-10245-MS, Los Alamos National Laboratory, Feb.
1985.
8. Ramos, J.I., Internal Combustion Engine Modeling, Newyork, 1989.
9. Abraham, J., Grasso, F., Wey, M.J and Bracco, F.V., Pressure Non-Uniformity
and Mixing Characteristics in Stratified-Charge Rotary Engine Combustion, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 880624, 1988.
10. Steinthorsson, E., Shih, T.I-P, Schock, H.j., and Stegman, J., Calculations of the
Unsteady, Three-Dimensional Flow Field Inside a Motored Wankel Engine, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 880625, 1988.
11. Abraham, J., and Bracco, F.V., Comparisons of Computed and Measured Mean
Velocity and Turbulence Intensity in a Motored Rotary Engine, SAE International
Congress and Exposition, Detroit, Mich., SAE Paper No. 881602, 1988.

70

12. Dimpelfeld, P. and Witze, P.O., Velocity Measurements in a 5.8l Stratified-Charge

Rotary Engine, Presented at the Fourth International Symposium on Applications of
Laser Anemometry to Fluid Mechanics, Lisbon, Portugal, July 11- 14, 1988.
13. Dimpelfeld, P.M., and Witze, P.O., Velocity Measurements in a 5.8l StratifiedCharge Rotary Engine: Tabulated Data, Sandia National Laboratories Report SAND888811, 1988.
14. Abraham, J., Bracco, F.V., Fuel-Air Mixing and Distribution in a Direct-Injection
Stratified-Charge Rotary Engine, SAE International Congress and Exposition, Detroit,
Mich., SAE Paper No. 890329, 1989.
15. Raju, M.S., and Willis, E.A., Analysis of Rotary Engine Combustion Process
Based on Unsteady, Three-Dimensional Computations, NASA Technical Memorandum
102469, AIAA-90-0643, 1990.
16. Hamady, F., Schock, H., and Stueckon, T., Air flow Visualization and LDV
Measurements in a Motored Rotary Engine Assembly-Part 1: Flow Visualization, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 900030, 1989.
17. Abraham, J., and Bracco, F.V., 3-D Computations of Premixed-Charge Natural
Gas Combustion in Rotary Engines, SAE International Congress and Exposition,
Detroit, Mich., SAE Paper No. 910625, 1991.
18. Abraham, J., and Magi, V., Effects of Ignition Cavity Flows on the Performance
of a Stratified-Charge Rotary Engine: Initial 3-D Predictions, SAE International
Congress and Exposition, Detroit, Mich., SAE Paper No. 941027, 1994.
19. Abraham, J and Bracco, F.V., 3-D Computations to Improved Combustion in a
Stratified-Charge Rotary Engine Part II: A Better Spray for the Pilot Injector, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 892057, 1989.
20. Abraham, J and Bracco, F.V., 3-D Computations to Improved Combustion in a
Stratified-Charge Rotary Engine Part III: Improved Ignition Strategies, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 920304, 1992.
21. Abraham, J., Bracco, F.V., and Epstein, P., 3-D Computations of Improved
Combustion in a Stratified-Charge Rotary Engine Part IV: Modified Geometries, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 930679, 1993.

71

22. Chouinard, E., Hamady, F., and Schock, H., Air flow Visulization and LDV
Measurement in a Motored Rotary Engine Assembly Part II: LDV Measurements, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 900031, 1989.
23. Abraham, J and Bracco, F.V., Comparison of Computed and Measured Pressure
in a Premixed-charge Natural-Gas-Fueled Rotary Engine, SAE International Congress
and Exposition, Detroit, Mich., SAE Paper No. 890671, 1989.
24. Roberts, J.A., Norman, T.J., Ekchian, J.A., and Heywood, J.B., Computer
Models for Evaluating Premixed and Disc Wankel Engine Performance, SAE
International Congress and Exposition, Detroit, Mich., SAE Paper No. 860613, 1986.
25. Shih, T, I-P., Schock, H.J., Ramos, J.I., Fuel-Air Mixing and Combustion in a
Two-Dimensional Wankel Engine, SAE International Congress and Exposition, Detroit,
Mich., SAE Paper No. 870408, 1987.
26. FLUENT Inc., FLUENT 6.1 documentation.
27. Sierens, R., Baert, R., Winterbone, D.E., and Baruah, P. C., A Comprehensive
study of Wankel Engine Performance, SAE International Congress and Exposition,
Detroit, Mich., SAE Paper No.830332, 1983.

72

APPENDIX
APPENDIX A - M-file used in the MATLAB for rotor and housing
M-file for rotor

% Rotor
R=0.1064;
%Generating Radius
e=0.01542;
%Eccentricity
C=0;
%Clearance
% v=(pi/6):.01:(pi/2);
% for the rotor face 1
% v=(5*pi/6):.01:(7*pi/6);
% for the rotor face 2
v=(3*pi/2):.01:(11*pi/6);
% for the rotor face 3
x=R*cos(2*v)-3*((e^2)/R)*sin(6*v).*sin(2*v)+2*e*(19*(e/R)^2*(sin(3*v)).^2).^(1/2).*cos(3*v).*cos(2*v);
y=R*sin(2*v)+3*((e^2)/R)*sin(6*v).*cos(2*v)+2*e*(19*(e/R)^2*(sin(3*v)).^2).^(1/2).*cos(3*v).*sin(2*v);
figure(1)
plot(x,y)
axis equal
hold on
M-file for housing
Zero degree rotation

% Rotor Housing
R=0.1064;
E=0.01542;
C=0.0015;
A=-pi/3:.01:pi/3;
t1=(0.01542)*cos(0*pi/6);
t2=(0.01542)*sin(0*pi/6);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+0*pi/6;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X+t1;
n=Y+t2;
figure(1)
plot(m,n)
axis equal
hold on

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

73

30 degree rotation

% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=-pi/2:.01:pi/6;
t1=(0.01542)*cos(1*pi/6-pi/2);
t2=(0.01542)*sin(1*pi/6-pi/2);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+1*pi/6;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X+t1;
n=Y+t2;
figure(1)
plot(m,n)
axis equal
hold on

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

60 degree rotation

Part A
% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=-pi/2:.01:0;
t1=(0.01542)*cos(2*pi/6);
t2=(0.01542)*sin(2*pi/6);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+2*pi/6;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X-t1;
n=Y-t2;
figure(2)
plot(m,n)
axis equal
hold on
Part B
% Rotor Housing
R=0.1064;

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

%Generating Radius

74

E=0.01542;
C=0.001;
A=pi/3:.01:pi/2;
t1=(0.01542)*cos(2*pi/6);
t2=(0.01542)*sin(2*pi/6);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+2*pi/6+pi;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X-t1;
n=Y-t2;
figure(2)
plot(m,n)
axis equal
hold on

%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

90 degree rotation

part A
% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=pi/6:.01:pi/2;
t1=(0.01542)*cos(3*pi/6-pi/2);
t2=(0.01542)*sin(3*pi/6-pi/2);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+3*pi/6+pi;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X-t1;
n=Y-t2;
figure(2)
plot(m,n)
axis equal
hold on
Part B
% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=-pi/2:.01:-pi/6;
t1=(0.01542)*cos(3*pi/6-pi/2);

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation

75

t2=(0.01542)*sin(3*pi/6-pi/2);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+3*pi/6;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X-t1;
n=Y-t2;
figure(2)
plot(m,n)
axis equal
hold on

%Eccentric rotation

120 degree rotation

part A
% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=-pi/2:.01:-pi/3;
t1=(0.01542)*cos(4*pi/6);
t2=(0.01542)*sin(4*pi/6);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+4*pi/6;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X+t1;
n=Y+t2;
figure(2)
plot(m,n)
axis equal
hold on
Part B
% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=0:.01:pi/2;
t1=(0.01542)*cos(4*pi/6);
t2=(0.01542)*sin(4*pi/6);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

76

theta=(atan(y./x))+4*pi/6+pi;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X+t1;
n=Y+t2;
figure(1)
plot(m,n)
axis equal
hold on
150 degree rotation

% Rotor Housing
R=0.1064;
E=0.01542;
C=0.001;
A=-pi/6:.01:pi/2;
t1=(0.01542)*cos(5*pi/6-pi/2);
t2=(0.01542)*sin(5*pi/6-pi/2);
x=E*cos(3*A)+(R+C)*cos(A);
y=E*sin(3*A)+(R+C)*sin(A);
r=sqrt(x.^2+y.^2);
theta=(atan(y./x))+5*pi/6+pi;
X=r.*cos(theta);
Y=r.*sin(theta);
m=X+t1;
n=Y+t2;
figure(2)
plot(m,n)
axis equal
hold on

%Generating Radius
%Eccentricity
%Clearance
%Eccentric rotation
%Eccentric rotation

In some cases two m-files like part A and B are used for the rotation of the housing.
This is due to the trigonometric limitations. For the remaining rotating angle it is just
the repeat of the above sequence.

77

APPENDIX B - Triangular valve

#include "udf.h"
#define NCOORDS 126 /* number of coordinates points*/
DEFINE_GEOM(slippery_arcs, domain, dt, position)
{
int i;
real dist2;
int closest_i;
real distmin =100.00;
real nodex[NCOORDS] = {
8.6603, 8.8987, 9.1151, 9.3087, 9.4793, 9.6262, 9.7493, 9.8481,
9.9224, 9.9720, 9.9969, 9.9969, 9.9720, 9.9224, 9.8481, 9.7493,
9.6262, 9.4793, 9.3087, 9.1151, 8.8987, 8.6603, 8.4003, 8.1194,
7.8183, 7.4978, 7.1587, 6.8017, 6.4279, 6.0380, 5.6332, 5.2144,
4.7825, 4.3388, 3.8843, 3.4202, 2.9476, 2.4676, 1.9815, 1.4904,
0.9957, 0.4985,-0.0000,-0.4985,-0.9957,-1.4904,-1.9815,-2.4676,
-2.9476,-3.4202,-3.8843,-4.3388,-4.7825,-5.2144,-5.6332,-6.0380,
-6.4279,-6.8017,-7.1587,-7.4978,-7.8183,-8.1194,-8.4003,-8.6603,
-8.8987,-9.1151,-9.3087,-9.4793,-9.6262,-9.7493,-9.8481,-9.9224,
-9.9720,-9.9969,-9.9969,-9.9720,-9.9224,-9.8481,-9.7493,-9.6262,
-9.4793,-9.3087,-9.1151,-8.8987,-8.6603,-8.4003,-8.1194,-7.8183,
-7.4978,-7.1587,-6.8017,-6.4279,-6.0380,-5.6332,-5.2144,-4.7825,
-4.3388,-3.8843,-3.4202,-2.9476,-2.4676,-1.9815,-1.4904,-0.9957,
-0.4985,-0.0000, 0.4985, 0.9957, 1.4904, 1.9815, 2.4676, 2.9476,
3.4202, 3.8843, 4.3388, 4.7825, 5.2144, 5.6332, 6.0380, 6.4279,
6.8017, 7.1587, 7.4978, 7.8183, 8.1194, 8.4003
};
real nodey[NCOORDS] = {
-5.0000,-4.5621,-4.1129,-3.6534,-3.1849,-2.7084,-2.2252,-1.7365,
-1.2434,-0.7473,-0.2493, 0.2493, 0.7473, 1.2434, 1.7365, 2.2252,
2.7084, 3.1849, 3.6534, 4.1129, 4.5621, 5.0000, 5.4255, 5.8374,
6.2349, 6.6169, 6.9824, 7.3305, 7.6604, 7.9713, 8.2624, 8.5329,
8.7822, 9.0097, 9.2148, 9.3969, 9.5557, 9.6908, 9.8017, 9.8883,
9.9503, 9.9876,10.0000, 9.9876, 9.9503, 9.8883, 9.8017, 9.6908,
9.5557, 9.3969, 9.2148, 9.0097, 8.7822, 8.5329, 8.2624, 7.9713,
7.6604, 7.3305, 6.9824, 6.6169, 6.2349, 5.8374, 5.4255, 5.0000,
4.5621, 4.1129, 3.6534, 3.1849, 2.7084, 2.2252, 1.7365, 1.2434,
0.7473, 0.2493,-0.2493,-0.7473,-1.2434,-1.7365,-2.2252,-2.7084,
-3.1849,-3.6534,-4.1129,-4.5621,-5.0000,-5.4255,-5.8374,-6.2349,
-6.6169,-6.9824,-7.3305,-7.6604,-7.9713,-8.2624,-8.5329,-8.7822,
-9.0097,-9.2148,-9.3969,-9.5557,-9.6908,-9.8017,-9.8883,-9.9503,
-9.9876,-10.0000,-9.9876,-9.9503,-9.8883,-9.8017,-9.6908,-9.5557,
-9.3969,-9.2148,-9.0097,-8.7822,-8.5329,-8.2624,-7.9713,-7.6604,
-7.3305,-6.9824,-6.6169,-6.2349,-5.8374,-5.4255
};
for(i=0; i < NCOORDS; i++)
{
dist2 = (position[0] - nodex[i])*(position[0] - nodex[i])
+ (position[1] - nodey[i])*(position[1] - nodey[i]);
if (dist2 < distmin)
{

78

closest_i = i;
distmin = dist2;

}
}
position[0] = nodex[closest_i];
position[1] = nodey[closest_i];

DEFINE_CG_MOTION(rotate, dt, vel,omega, time, dtime)

{
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);
omega[2] = 1.0;
}

79

APPENDIX C UDF for rotary engine

#include "udf.h"
#define NCOORDS 290 /* number of coordinates points*/
#define R -0.01542
DEFINE_GEOM(slippery_arcs, domain, dt,position)
{
int i;
real dist2;
int closest_i;
real distmin = 100.00;
real nodex[NCOORDS] = {
-0.1258,-0.1258,-0.1256,-0.1253,-0.1249,-0.1244,-0.1237,-0.1230,
-0.1221,-0.1212,-0.1201,-0.1189,-0.1177,-0.1163,-0.1148,-0.1133,
-0.1116,-0.1099,-0.1081,-0.1063,-0.1043,-0.1023,-0.1003,-0.0981,
-0.0960,-0.0938,-0.0915,-0.0892,-0.0869,-0.0846,-0.0822,-0.0798,
-0.0774,-0.0750,-0.0726,-0.0702,-0.0678,-0.0654,-0.0630,-0.0606,
-0.0583,-0.0559,-0.0536,-0.0514,-0.0491,-0.0469,-0.0447,-0.0426,
-0.0405,-0.0384,-0.0364,-0.0344,-0.0325,-0.0306,-0.0287,-0.0269,
-0.0251,-0.0233,-0.0216,-0.0200,-0.0183,-0.0167,-0.0152,-0.0136,
-0.0121,-0.0106,-0.0092,-0.0077,-0.0063,-0.0049,-0.0035,-0.0021,
-0.0007, 0.0007, 0.0021, 0.0035, 0.0049, 0.0063, 0.0077, 0.0092,
0.0106, 0.0121, 0.0136, 0.0152, 0.0167, 0.0183, 0.0200, 0.0216,
0.0233, 0.0251, 0.0269, 0.0287, 0.0306, 0.0325, 0.0344, 0.0364,
0.0384, 0.0405, 0.0426, 0.0447, 0.0469, 0.0491, 0.0514, 0.0536,
0.0559, 0.0583, 0.0606, 0.0630, 0.0654, 0.0678, 0.0702, 0.0726,
0.0750, 0.0774, 0.0798, 0.0822, 0.0846, 0.0869, 0.0892, 0.0915,
0.0938, 0.0960, 0.0981, 0.1003, 0.1023, 0.1043, 0.1063, 0.1081,
0.1099, 0.1116, 0.1133, 0.1148, 0.1163, 0.1177, 0.1189, 0.1201,
0.1212, 0.1221, 0.1230, 0.1237, 0.1244, 0.1249, 0.1253, 0.1256,
0.1258, 0.1258, 0.1258, 0.1256, 0.1253, 0.1249, 0.1244, 0.1237,
0.1230, 0.1221, 0.1212, 0.1201, 0.1189, 0.1177, 0.1163, 0.1148,
0.1133, 0.1116, 0.1099, 0.1081, 0.1063, 0.1043, 0.1023, 0.1003,
0.0981, 0.0960, 0.0938, 0.0915, 0.0892, 0.0869, 0.0846, 0.0822,
0.0798, 0.0774, 0.0750, 0.0726, 0.0702, 0.0678, 0.0654, 0.0630,
0.0606, 0.0583, 0.0559, 0.0536, 0.0514, 0.0491, 0.0469, 0.0447,
0.0426, 0.0405, 0.0384, 0.0364, 0.0344, 0.0325, 0.0306, 0.0287,
0.0269, 0.0251, 0.0233, 0.0216, 0.0200, 0.0183, 0.0167, 0.0152,
0.0136, 0.0121, 0.0106, 0.0092, 0.0077, 0.0063, 0.0049, 0.0035,
0.0021, 0.0007,-0.0007,-0.0021,-0.0035,-0.0049,-0.0063,-0.0077,
-0.0092,-0.0106,-0.0121,-0.0136,-0.0152,-0.0167,-0.0183,-0.0200,
-0.0216,-0.0233,-0.0251,-0.0269,-0.0287,-0.0306,-0.0325,-0.0344,
-0.0364,-0.0384,-0.0405,-0.0426,-0.0447,-0.0469,-0.0491,-0.0514,
-0.0536,-0.0559,-0.0583,-0.0606,-0.0630,-0.0654,-0.0678,-0.0702,
-0.0726,-0.0750,-0.0774,-0.0798,-0.0822,-0.0846,-0.0869,-0.0892,
-0.0915,-0.0938,-0.0960,-0.0981,-0.1003,-0.1023,-0.1043,-0.1063,
-0.1081,-0.1099,-0.1116,-0.1133,-0.1148,-0.1163,-0.1177,-0.1189,
-0.1201,-0.1212,-0.1221,-0.1230,-0.1237,-0.1244,-0.1249,-0.1253,
-0.1256,-0.1258
};
real nodey[NCOORDS] = {
-0.0000,-0.0034,-0.0068,-0.0102,-0.0135,-0.0169,-0.0202,-0.0235,

80

-0.0267,-0.0299,-0.0331,-0.0362,-0.0392,-0.0422,-0.0452,-0.0480,
-0.0508,-0.0535,-0.0562,-0.0587,-0.0612,-0.0636,-0.0659,-0.0681,
-0.0703,-0.0723,-0.0743,-0.0761,-0.0779,-0.0796,-0.0811,-0.0826,
-0.0840,-0.0853,-0.0865,-0.0877,-0.0887,-0.0897,-0.0906,-0.0914,
-0.0921,-0.0927,-0.0933,-0.0938,-0.0943,-0.0947,-0.0950,-0.0953,
-0.0955,-0.0957,-0.0959,-0.0960,-0.0960,-0.0961,-0.0961,-0.0961,
-0.0961,-0.0960,-0.0959,-0.0959,-0.0958,-0.0957,-0.0956,-0.0955,
-0.0954,-0.0953,-0.0952,-0.0952,-0.0951,-0.0951,-0.0950,-0.0950,
-0.0950,-0.0950,-0.0950,-0.0950,-0.0951,-0.0951,-0.0952,-0.0952,
-0.0953,-0.0954,-0.0955,-0.0956,-0.0957,-0.0958,-0.0959,-0.0959,
-0.0960,-0.0961,-0.0961,-0.0961,-0.0961,-0.0960,-0.0960,-0.0959,
-0.0957,-0.0955,-0.0953,-0.0950,-0.0947,-0.0943,-0.0938,-0.0933,
-0.0927,-0.0921,-0.0914,-0.0906,-0.0897,-0.0887,-0.0877,-0.0865,
-0.0853,-0.0840,-0.0826,-0.0811,-0.0796,-0.0779,-0.0761,-0.0743,
-0.0723,-0.0703,-0.0681,-0.0659,-0.0636,-0.0612,-0.0587,-0.0562,
-0.0535,-0.0508,-0.0480,-0.0452,-0.0422,-0.0392,-0.0362,-0.0331,
-0.0299,-0.0267,-0.0235,-0.0202,-0.0169,-0.0135,-0.0102,-0.0068,
-0.0034, 0.0000, 0.0034, 0.0068, 0.0102, 0.0135, 0.0169, 0.0202,
0.0235, 0.0267, 0.0299, 0.0331, 0.0362, 0.0392, 0.0422, 0.0452,
0.0480, 0.0508, 0.0535, 0.0562, 0.0587, 0.0612, 0.0636, 0.0659,
0.0681, 0.0703, 0.0723, 0.0743, 0.0761, 0.0779, 0.0796, 0.0811,
0.0826, 0.0840, 0.0853, 0.0865, 0.0877, 0.0887, 0.0897, 0.0906,
0.0914, 0.0921, 0.0927, 0.0933, 0.0938, 0.0943, 0.0947, 0.0950,
0.0953, 0.0955, 0.0957, 0.0959, 0.0960, 0.0960, 0.0961, 0.0961,
0.0961, 0.0961, 0.0960, 0.0959, 0.0959, 0.0958, 0.0957, 0.0956,
0.0955, 0.0954, 0.0953, 0.0952, 0.0952, 0.0951, 0.0951, 0.0950,
0.0950, 0.0950, 0.0950, 0.0950, 0.0950, 0.0951, 0.0951, 0.0952,
0.0952, 0.0953, 0.0954, 0.0955, 0.0956, 0.0957, 0.0958, 0.0959,
0.0959, 0.0960, 0.0961, 0.0961, 0.0961, 0.0961, 0.0960, 0.0960,
0.0959, 0.0957, 0.0955, 0.0953, 0.0950, 0.0947, 0.0943, 0.0938,
0.0933, 0.0927, 0.0921, 0.0914, 0.0906, 0.0897, 0.0887, 0.0877,
0.0865, 0.0853, 0.0840, 0.0826, 0.0811, 0.0796, 0.0779, 0.0761,
0.0743, 0.0723, 0.0703, 0.0681, 0.0659, 0.0636, 0.0612, 0.0587,
0.0562, 0.0535, 0.0508, 0.0480, 0.0452, 0.0422, 0.0392, 0.0362,
0.0331, 0.0299, 0.0267, 0.0235, 0.0202, 0.0169, 0.0135, 0.0102,
0.0068, 0.0034
};
for(i=0; i < NCOORDS; i++)
{
dist2 = (position[0] - nodex[i])*(position[0] - nodex[i])
+ (position[1] - nodey[i])*(position[1] - nodey[i]);
if (dist2 < distmin)
{
closest_i = i;
distmin = dist2;
}

}
position[0] = nodex[closest_i];
position[1] = nodey[closest_i];

DEFINE_CG_MOTION(eccentric, dt, vel, omega, time, dtime)

{
real beta;
NV_S (vel, =, 0.0);
NV_S (omega, =, 0.0);

81

omega[2] = 1.0;
beta = 3;
vel[0] = -beta*R*sin(beta*time);
vel[1] = beta*R*cos(beta*time);
vel[2] = 0;
}

82