SCHOOL OF MECHANICAL AND

MANUFACTURING ENGINEERING

FINITE ELEMENT METHOD

Final Project
Analysis and Optimization
of Energy Storage Flywheel with FEM
by

Minh Nguyen Vo Nhat

Student ID:

z5026422

Faculty of Engineering
April 2015

EXECUTIVE SUMMARY
1/ Problem
This report will use FEM to analyze three typical flywheel used for flywheel energy storage system
which the shape is hollow solid cylinder. The dimension of the flywheel is constant 0.1 x 0.2 x 0.06 (inner
radius x outer radius x height).

Objective: analyse and conclude a optimised choice of material and structure for a flywheel
used for flywheel energy storage system

2/Case 1: Isotropic material:

Below is the model the the mesh of the first case

We have the stresses along radius and equivalent stresscontour plot:

After calculate, we know maximum rotational velocity and maximum energy can be stored:
"_$%&'(&)$* =>

1
Ek = h (ro 4 ri 4 ) 2 =37954.61 J
4
2

2/Case 2: Isotropic material:

st

This case use the same model the the mesh as 1 case, the diferent is materiAL
We have the stresses along radius and equivalent stresscontour plot:

After calculate, we know maximum rotational velocity and maximum energy can be stored:
"_$%&'(&)$* "_+_,-./(_&('0&'(&)$* =47091.2 J
2/Case 1: Isotropic material:
Below is the model the the mesh of the first case

We have the stresses along radius and equivalent stresscontour plot:

After calculate, we know maximum rotational velocity and maximum energy can be
stored:
"_$%&'(&)$* =>

1
Ek = h (ro 4 ri 4 ) 2 =37954.61 J
4

This report includes multiple validations of geometry of numerical model: meshing

quality convergence, 3D and 2D comparison; and stresses results verifications: by theory for

one rim isotropic and orthotropic case, by published paper for multi-rims flywheel case.
Below is stress results compared with a published paper:

After calculation of 3 cases, we have a table to sum it up:

One-rim

Case

Aluminum Alloy

Maximum velocity can withstand

(rad/s)

1453.25

One-rim T300

Two-rim T300

2066.096

5564.14

Maximum energy capacity (J)

37954J

47091.2

341549.27

Mass = (Kg)

5.155524

3.294324

3.294324

Table 1. Results of 3 analyzed cases

Although there was a limitation, the model was validated with multi approach for each
of case. The report is able to conclude that multi-rim composite flywheel have the best
performance as predicted. So that is the optimized choice for materials and structure of a
flywheel using for energy storage systems.
In general, all of two objectives is achieved. Due to time and size limit of this paper,
the variation of each case was not examined, like changing some material, dimension
.
These will be next objectives in a future report. Thank you for reading.

CONTENTS
INTRODUCTION

Background

Problem Description

Assumption

Objective

METHODOLOGY

Theory and formulas

2.1.1.

Kinetic energy of flywheel

2.1.2.

Stress analysis of flywheel, isotropic case:

10

2.1.3.

Stress analysis of flywheel, orthotropic case:

10

Single-rim flywheel by isotropic material (Aluminium)

11

Single-rim flywheel with orthotropic elasticity composite material

(T300s/Epoxy),

14

Two-rim flywheel with orthotropic elasticity composite material

(T300s/Epoxy),

17

VALIDATION

20

Verification of numerical model geometry

20

3.1.1.

3D or 2D?

20

3.1.2.

Mesh convergence study

21

Stress results validation

21

3.2.1.

Isotropic case:

21

3.2.2.

Single-rim orthotropic case:

22

Multi-rim orthotropic case:

22

LIMITATIONS 24
CONCLUSION 25
REFERENCE

26

APPENDIX A: Orthotropic stresses calculation code MATLAB of Kamf (2012)

27

INTRODUCTION

Background
Energy storage assumes is a huge issue to todays society, where nearly all that
we use in our life need energy. If the energy can't be directly supplied to where it is
needed, it may be taken from some sort of storages. For example, it is wood,
gasoline, or different kind of batteries.
Flywheel is one of energy storage type, which was employed in the ancient time.
It has been developed to many generations and type. One of them is composite
flywheel. In the late 1970s and early 1980s, much research on composite flywheels
was carried out in response to the energy storage needs in space and automotive
applications, Many practical energy storage devices which use composite fly-wheels,
for example, non-intermittency power supplies (Andoh et al. 1996), batteries for
automobiles (Jost 1994), and pulsed power supplies for electromagnetic guns
(Curtiss et al. 1995). For these flywheel energy storage systems to be competitive
with electrochemical batteries, it is essential that they contain high energy density
flywheels. The below figure is from KERS (Flybrid System 2010), which is a flywheel
energy storage system can store maximum 60 KWh use for F1 car.

Figure 1. KERS (Flybrid System 2010)

Then Chamis, C and Kiraly L (1976) from NASA and numbers of research have
conduct an analysis about two type of composite flywheel: single-rim and multi-rim.
This is for utilizing the radical stress of flywheel, and the rims are assembled by
interference fit.

Problem Description
This report will use FEM to analyze three typical flywheel used for flywheel
energy storage system which the shape is hollow solid cylinder. The dimension of the
flywheel is constant 0.1 x 0.2 x 0.06 (inner radius x outer radius x height).

Single-rim flywheel by isotropic material (Aluminium), Figure 2.a

Single-rim

flywheel

with

orthotropic

elasticity

composite

material

elasticity

composite

material

(T300s/Epoxy), Figure 2.b

Two-rim

flywheel

with

orthotropic

(T300s/Epoxy), Figure 2.c

(a)

(b)

(d)

(c)

Figure 2. Flywheel; (a) Single-rim flywheel by isotropic material flywheel, (b) Single-rim
flywheel with orthotropic elasticity composite material, (c) Two-rim flywheel with
orthotropic elasticity composite material, (d) Multi-rim flywheel model( Ha et al. 2012)
Parameters

Aluminium

T300/Epoxy

Longitudinal Youngs modulus

Radial Youngs modulus

Er

Longitudinal tensile strength

(MPa)

(MPa)

(MPa)

Longitudinal compressive strength

'

(MPa)
Radial tensile strength

Radial compressive strength

Poisson ratio

(MPa)

r (MPa)

r / z / rz

r / z / rz
Density (kg/m 3 )

Shear strength

(Mpa)
(MPa)

7.1E5

1.81E5

7.1E5

1.03E5

280

1500

280

1500

280

40

280

246

0.33/0.33/0.33

0.3/0.25/0.25

2E5

2/2/5 E3

2770

1600

Table 2. Properties of materials

Assumption
The following assumptions were applied to keep the hand calculation more
simple, this is in response for more analysis by Finite Element Method:

In all equations an assumption of a system in thermal equilibrium is

made, meaning that no temperature change over time is present.

A rigid installation is used to connect the flywheel and its drive-shaft, no

key-ways are needed to fix/drive the flywheel and there exist no slide,
built-in stress and deformations on the connection surface; therefore,
displacement constraints can be simply applied on the shaft hole;

The flywheel only works in the vertical plane (X-Y plane) so that the
gravity can be simply applied;

The aero-dynamical resistance can be neglected;

There exists no vibration;

The fillets/chamfers can be neglected unless dimensioned;

Objective

These following objectives are the expected outcome of this project, and the
scope of this project is not limited:

Analyse three typical type of flywheel base on maximum energy before failure

Optimise the flywheel by concluding the suitable material and structure for high
energy storage flywheel base of energy density.
So these upper goals are very important for manufacturing flywheel using for

energy storage and have big research potential.

METHODOLOGY
This project analyses 3 different cases but to understand and validate, there
are two sections: theory and FEA

Theory and formulas

2.1.1. Kinetic energy of flywheel
The total energy stored in a rotating mass can be calculated by using:

1
Ek = I 2
2
With

(2.1)

I is inertia moment (Kg.m2)

is rotational velocity ( rad/s)

The moment of inertia of a barrel cylinder can then be calculated using :

I=

1
1
m ( ri 2 + ro 2 ) = V ( ri 2 + ro 2 )
2
2

1
h (ro 2 ri 2 )( ri 2 + ro 2 )
2
1
= h (ro 4 ri 4 )
2
=

(2.2)

And substitute (1.2) to (1.1), we have:

1
=> Ek = h (ro 4 ri 4 ) 2
4

(2.3)

So the energy is proportional with density and rotational velocity. But the
problem here is how fast the flywheel can withstand. This is up to the stresses due to
the rotation

2.1.2. Stress analysis of flywheel, isotropic case:

The general equation of motion is:

. + F = u

(2.4)

For a cylinder coordinate, stress is described as the figure below:

Figure 3. Cylinder coordinate description and differential piece definition (Kamf 2012)
As a classical problem, Hearn (1997) has derived (2.5) for the case of isotropic
materials:

ri 2ro 2
2
r (r ) = (3 + v)
r
+
r
+

r
i
o

8
r2

(2.5)

And

2
ri 2ro 2
2
2
(r ) =
(3
+
v
)
r
+
r
+

(1
+
3
v
)
r

o
8
r2

(2.6)

2.1.3. Stress analysis of flywheel, orthotropic case:

The equations (2.6) and (2.7) are enough for most basic case of isotropic
flywheel. But for composite, which is assumed orthotropic, is much more complex. As
Kamf (2012) have summarize previous work. These are equations for stress
calculation:

r (r ) = C1r

(1+ )

+ C1r

(1 )

10

2r 2
+ 2
(3 + )
9

(2.7)

 (r ) = (C1r

(1 )

+ C1r

(1 )

2r 2
) + 2 9 (3 + )

(2.8)

With

C1 =

(rr
K ( ri 3ro ro 3ri )
i o)

) ; C = ( K (r
2

(ri 2 ro 2 )

+3

ri +3 )

(ri 2 ro 2 )

E 1 vzr vrz
v r + vzr v z
2
;K =
=
; =
(3 + )
Er 1 vz v z
1 vz v z
2 9
2

Single-rim flywheel by isotropic material (Aluminium)

First, a FEA program named ANSYS was used to calculate the stresses of
rotating cylinder. The FE model is on 2D with axisymmetric by Y behavior, the cylinder
was model as a haft of cross section as a rectangle, which will be rotating with given
dimension and input as table below.
Dimension (m)

0.1 x 0.2 x 0.06

Rotation velocity

5000 rad/s

Axisymmetric axis

Element types

PLANE183

Element size (m)

0.005

Table 3. 2D FE model properties

10 mm

60 mm

Rotational Velocity

Flywheel Cross Section

Axisymmetric Axis

100 mm

Figure 4. Problem model

And Figure 4 is 2D FE model with only boundary conditions: axisymmetric and rotational
speed. And the meshed model with Mapped Method as was stated in Figure 5:

11

Figure 5. 2D FE mesh model

After solving, plot nodal contour result of normal stress on X, Y and Z axis, we
have Figure 6. Because this is 2D model, so we have stresses convert to cylinder
coordination:

X = r ; Y = z ; X =

(2.9)

It is a easy way to consider failure in either isotropic, when the materials

properties is identical in every direction, or in cylindrical coordinator. The reason why
we must normalize the result in a cylindrical is: for each cross-section segment along
a full 3D cylinder disk, 3 of cylindrical stresses are equal to other cross-section
segments. Surprisingly, in a Beta Option of ANSYS, there was a Visual Expansion
which can express a full expand or quartered expand. This is just help for visual but
not have a real result of a 3D model. By this, I decided to build another model but on
3D to validate the result of 2D axisymmetric model in section 2.2.2.1. If only the
results are identical, this 2D model will be use thought out the report due to its
compactness of the 2D model and simplicity of mesh control.

12

Figure 6. 2D FEA von-Mises equivalent stress results when =5000 rad/s

We have result of von-Mises equivalent stress as the following Figure 6

According to the plot in Figure 7 below, the tangential stress 6 distributed in a
pattern of decreasing along the radius, as the point moving close to the inner radius,
the stress is maximum or vice versa. While radial stress ( distributed in a pattern of
a up-side-down hyperbolic, the maximum radial stress in near the middle of the cross
section. As the stress along the z axis 7 is so small, which can be negligible, and the
radial stress is also small, the Von-mises equivalent stress is basically up to
tangential stress. We can also see the pattern and amplitude of the Von-mises
equivalent stress is look like these of tangential stress.

Aluminium, omega=5000rad/s
800.00

Radial Stress
Tangential Stress
Axial Stress
Equivalent Stress

600.00
MPa

400.00
200.00
0.00
0

12.5

25

37.5

50

62.5

75

87.5

100

-200.00
mm

Figure 7. Stresses along radius of isotropic material flywheel when =5000

rad/s

Linearity of model

Hence, in this isotropic case, the factor most effects to energy capacity is
tangential stress, which was calculate by (2.7) and it is approximately simplified by:

6 = 6 :

(2.10)

By running simutaneously this FE model with different velocity. We can calculate

6 and prove this linear equation is right. And this lead to we can calculate the
maximum energy, which respect to

(Equation 2.3), this model can store.

13

We have a table as follow:

Velocity

Tangential Stress

Variation

(rad/s)

(MPa)

of A

5000

2.34

2.8296E-05

0.00%

6000

3.34

2.83665E-05

-0.08%

7000

4.48

2.84504E-05

-0.09%

8000

5.76

2.85478E-05

-0.10%

9000

7.17

2.86592E-05

-0.12%

10000

8.67

-0.13%
2.8785E-05
Table 4. Differents results with increasing the velocity of isotropic case
As the result above, we can see this reports linear approximation (2.11) is true
with aceptable error (<1%). By this, we have a calculation of maximum velocity that a
alluminium flywheel can withstand

;-< =

>?_@ABCD
E?

:FG
:.F:IJKLGM

= 1453.25

(-S
%

(2.11)

For this isotropic flywheel, from (2.7) maximum energy is:

"_$%&'(&)$* =>

1
Ek = h (ro 4 ri 4 ) 2 =37954.61 J
4

Single-rim flywheel with orthotropic elasticity composite material

(T300s/Epoxy),
As we concluded in the isotropic section, the energy capacity is depend on
tangential yield stress.

So to be suitable, when composite have more strength in

desired direction, we can choose a wound filament composite flywheel. With the
orientation of 900, the longitudinal direction of the fiber is now the tangential direction
of the flywheel.
Isotropic materials have the properties identical in every direction while
orthotropic materials have its properties different in a normal Descartes xyz
coordinator. However, the orthotropic have cylindrical properties. That is at every
point of the material, the properties are symmetric by 3 axis: axial, radial and
circumferential. In the stress result of axisymmetric is also have the same properties.

14

In general, we can compare the stresses of the axisymmetric model with cylindrical
properties of orthotropic material.
For the setup and FE model, it is the same as isotropic case, just change the
material to T300 and applying cylindrical normalization to it. We have the result of
Von-Mises equivalent stress as Figure 8:

Figure 8. 2D FEA von-Mises equivalent stress results of one-rim composite

flywheel
According to the plot in Figure 9, the tangential stress 6 and radial stress (
have the same pattern in the isotropic case. However, the axial stress along the z
axis s ( is now having noticeable amplitude which is pointed close to the inner
radius. But again, the axial stress and radial stress is also small, the von-Mises
equivalent stress is basically up to tangential stress. We can also see the pattern
and amplitude of the von-Mises equivalent stress is look like tangential stresses, as
the isotropic case above.
Single-rim T300, omega=5000rad/s
800.00
Radial Stress
Tangential Stress
Axial Stress
Equivalent Stress

600.00

MPa

400.00
200.00
0.00
0

12.5

-200.00

25

37.5

50

62.5

75

87.5

100

mm

Figure 9. Stresses along radius of one-rim orthotropic flywheel =5000 rad/s

Linearity of model

15

In this case, due to the strength is different on 3 normalized direction, we must

consider both the radial stress and tangential stress causes, which was calculate by
(2.7) and (2.8) and it is approximately simplified by:

6 = 6 :

(2.12)

( = ( :

(2.13)

Again, by running simutaneously this FE model with different velocity. We can

calculate 6 and ( and prove this linear equation (2.13) and (2.14) is right.
We have the table as below, which is stresses coresspondent to particular
velocity, and 2 value of A was calculated by interpret (2.13) and (2.14):
Velocity

Radial

Tangential

(rad/s)

Stress (MPa)

Stress (MPa)

5000

234.26

6000
7000
8000
9000
10000

715.79

2.86E-05

333.63

1030.22

448.04

Variation of

Variation of

9.37E-06

0.00%

0.00%

2.86E-05

9.37E-06

0.05%

0.01%

1401.43

2.86E-05

9.37E-06

0.11%

0.05%

575.89

1829.31

2.86E-05

9.36E-06

0.17%

0.12%

716.62

2313.74

2.86E-05

9.35E-06

0.23%

0.17%

867.43

2854.85

2.85E-05

9.35E-06

0.29%

0.21%

Table 5. Differents results with increasing the velocity of one-rim orthotropic case
As the result above, we can see this reports linear approximation (2.11) is true
with aceptable error of less than 1% when speed is smaller than 50000 rad/s(which
is very fast speed). By this, we have a calculation of maximum velocity that a
alluminium flywheel can withstand

;-< =
;-< =

>T_@ABCD
ET
>?_@ABCD
E?

+MGG
:.FJKLGM
YG
I.Z[KLGJ

= 7238.06

= 2066.096

(-S
%
(-S
%

(2.14)
(2.15)

As we can see , the maximum velocity is up to which is smaller in (2.16) and

(2.17). For this single-rim orthotropic flywheel, from (2.7) maximum energy is:
"_+_,-./(_&('0&'(&)$* =47091.2 J
The maximum velocity repect to tensile radial stresses is smaller than 3.6 times
when compare with tangential stress. The result is the energy is lower by 3.6^2=10.6
times. However, it is higher than isotropic case, which is just 37954J.

16

So we must find an approach to reduce radial stress of othortropic case. That

is layer-up this material in section 2.4

Two-rim flywheel with orthotropic elasticity composite material

(T300s/Epoxy),
As Ha et al (2009) claims, the efficiency of composite flywheel can be increase
by layer up to multi-layer, but in this report, only two-layered flywheel was analyzed.
The following figure 10 is the FE model of it:
Rotational

10

60 mm

Velocity

Axisymmetric Axis

Layer 1
100 mm

Layer 2
50 mm

Figure 10. Problem model of two-rim flywheel

This model uses 2 layer of T300. It was model as 2 separate flywheels, which
is each of the layers. However, when the flywheel is spinning, the layers displace and
expand their diameter. And now the flywheel is disengaged and the outer layers
doesnt have any torque to make they spin. Normally, they will apply interference fit to
each layer, which means pre-stressed (Ha et al, 2008). But because we looking to
only the case when it reaches maximum allowance speed, and then every layer prestressed status is nearly released or their contacts is nearly separate, and can be
model as separate layers as Figure 10. So we have an assumption for this case only
and this assumption will be discussed later:

Asumption: The rims have friction contact with other all the time.

After mesh with the same properties in Table 6, we have:

17

Figure 11. 2D FE mesh model

Contacts define: As assume that 2 layers is separate, a friction contact

type with 0.2 coefficient was used between two layer. All other option is
Program Controlled

One of the contour results as Figure 10:

Figure 12. Tangential stress of two-rim flywheel =5000 rad/s

According to Figure 13 below, the tangential stress and radial stress divide into
2 parts as 2 parts of the flywheel model. As we can see, the tangential stress in
inner rim is smaller than one of outer rim. The maximum tangential stress now is just
less than half of one layer case. And the radial stresses of 2 rims is roughly similar,
however, it is ten times smaller when compare with one layer case. This mean the
flywheel is able to spin faster without reaching yield stress. So with this layer up has
solved the problem of one layer case.

18

Double-rim T300, omega=5000rad/s

300.00

Tangential Stress
Axial Stress
Radial Stress

MPa

200.00
100.00
0.00
0

12.5 25 37.5 50 62.5 75 87.5 100

-100.00

mm

Figure 13. Stresses along radius of two-rim orthotropic flywheel =5000 rad/s

Linearity of model:

As the result of one-layer orthotropic, we can see this system would have linear
behavior. The reason is just imagine this two-layer flywheel is 2 separate one-layer
flywheels. Again, by running simutaneously this FE model with different velocity. We
can calculate 6 and ( and prove linear equations (2.13) and (2.14) is right.
Velocity

Radial

Tangential

Stress

Stress

(rad/s)

(MPa)

(MPa)

5000

32.30

257.24

1.292E-06

6000

46.48

370.38

7000

63.25

8000

Variation

Variation

of

of

1.029E-05

0.00%

0.00%

1.291E-06

1.028E-05

0.09%

0.01%

503.95

1.291E-06

1.028E-05

0.12%

0.05%

82.58

658.03

1.290E-06

1.028E-05

0.12%

0.12%

9000

104.45

832.10

1.289E-06

1.027E-05

0.19%

0.17%

10000

128.93

1027.20

1.288E-06

1.027E-05

0.26%

0.21%

Table 7. Differents results with increasing the velocity of one-layer orthotropic case
By this, we have a calculation of maximum velocity that a two-layer flywheels
flywheel can withstand
^_` =
^_` =

ab_cdefg
hb
al_cdefg
hl

=
=

+MGG
+.G:FKLGM
YG
+.:I:mLGJ

19

= 12073.63
= 5564.14

i_j
k
i_j
k

(2.16)
(2.17)

As we can see, the maximum velocity is 5564.14

(-S
%

. For this two-rim orthotropic

flywheel, from (2.7) maximum energy is:

"_:_,-./(_&('0&'(&)$* =341549.27 J=341 KJ
As we can see with the result, layer up 2-rim of material lowered the stresses,
both the radial and tangential one. So this case is the best among 3 cases which
was analyzed.

VALIDATION
Verification of numerical model geometry
3.1.1. 3D or 2D?
Due to validate with more real simulation, a 3D model was analyzed to make
sure the result is the same, so that the 2D model can be used for further study. The
3D model is 1/8 of a full flywheel, all properties is the same as Aluminum case
(Section 2.2) with same velocity, just different about element type and axisymmetric
axis. And we use this command to get ANSYS know this is 1/8 of full model:
/prep7
cyclic,8
/solu
Dimension (m) (outer radius/ inner radius/ height)

0.2/0.1/0.06

Rotation velocity

5000 rad/s

Axisymmetric axis

Element type

SOLID186

Element size (m)

0.005

Table 8. 3D FE model properties

And this is 3D FE model and it have exactly the same result as 2D:

20

Figure 14. 3D FE model of an eighth of a full flywheel left ; 3D FEA stress results of
tangential stress middle; full result expandsion - right

So we can use 2D model for a faster calculation and an extreme convergence

mesh quality as presented in following section.

3.1.2. Mesh convergence study

In order to check the quality of current mesh, a convergence study of mesh
sizing was conducted. This study uses the aluminum flywheel to calculate, with the
velocity of 5000 rad/s. In the table below, the smaller element size or the mesh is
finer, but the result doesnt change. The Mapped Face Meshing control turns the
geometry to smaller rectangles. And structure to totally the same, this lead to this
extremely convergence mesh.
Edge

Number

Stress

Sizing

of

(mm)

Element

2591

707.87042236328125

0%

2640

707.87042236328125

0%

2795

707.87042236328125

0%

0.5

3084

707.87042236328125

0%

0.2

4020

707.87042236328125

0%

0.1

5824

707.87042236328125

0%

Maximum Stress(Pa)

Variation
ratio (%)

Table 9. Convergence study of mesh quality

Stress results verification

3.2.1. Isotropic case:
To clearly see the stress distribution along radius, a plot was made to compare
FEA and theoretical result from equations (2.6) and (2.7). It is quite a simple formula
to calculate by Excel:

21

800.00

Aluminium, omega=5000rad/s
FEA Radial Stress

600.00

FEA Tangential Stress

400.00

200.00

0.00
0

12.5

25

37.5

50

62.5

75

87.5

100

Figure 15. Theoretical result and FEM result isotropic case =5000 rad/s
As we can see, the result in theory is similar to FEA, however, it is bigger around
1% because of theoretical formulas (2.6) and (2.7) assume axial stresses is zero. So
the result is acceptable.

3.2.2. Single-rim orthotropic case:

By the MATLAB code of Kamf (2012) base on (2.8) and (2.9) in Appendix
B, we can calculate the stresses of this orthotropic case, the results was plotted as
below:
800.00

One-rim T300, omega=5000rad/s

FEA Radial Stress
FEA Tangential Stress
Theoretical Radial Stress

600.00
400.00
200.00
0.00
0

12.5

25

37.5

50

62.5

75

87.5

100

Figure 16. Theoretical result and FEM result orthotropic case =5000 rad/s
Again, the result in theory is similar to FEA and it is bigger around 1.5% because
of theoretical formulas (2.6) and (2.7) assume axial stress is zero. So the result is
acceptable.

Multi-rim orthotropic case:

22

To validate this case, the author uses a published paper by Wen, S and Jiang, S
in 2012. They did a research on multi-rim flywheel using a combination of material
including: aluminum, e-glass and T300/2500. This paper has interference fit effects
on result. The figure and table below are their flywheel model and the result of this
paper.
Rim
Radius(In-Out)

Inner rim

Middle rim

Outer rim

(Aluminum)

(E-Glass)

(T300)

0.06-0.08

0.06-0.08

0.06-0.08

Table 10. Wen, S and Jiang, S Multi-rim flywheel configuration

Figure 17. Analytical result when in maximum velocity = 58114 RPM

The author builds the same model on ANSYS, with the same method as Section
2.4. The result is roughly the similar in shape and amplitude.
90
Radial Stress

Axis Title

1800
1500
1200
900
600
300
0

MPa

60
30

Tangential Stress
60

82.5

105

127.5

150

0
60

mm

82.5

105
mm

23

127.5

150

Figure 17. FEA Tangential stress and radial stress result when = 58114 RPM top;
Tangential stress contour plot - bottom
So with the comparison with a published paper, the multi-rim flywheel FE model
is validated, and proves a not-considering- interference-fits model has a correct result
when it reaching maximum velocity.

LIMITATIONS
There was a limitation of this paper that the interference fit effect for
assembling 2 rims together wasnt conduct successfully so that the author cant take
these effects into account. Thats why assumption in section (2.4) was made, in fact,
there was a gap between rims or there wasnt any contact when it is spinning. The
author tried to use non-linear offset contacts to simulate the interference fit. It shows
a good result when there was no velocity. But when applied a rotational velocity, the
solutions is not convergence and appear errors as Figure 19.

Figure 17. Nonlinear contacts error when applied velocity

This is the result of interference fit with no rotaion, using nonlinear contact:

24

Figure 18. Equivalent stress when applied interference fit with no rotation
And the figure 18 below is the gap in the case of Section 2.4:

Figure 18. Gap in between 2 layer in 2.4 case.

According to various published experiments and calculations (Kamf 2012, Ha et
al 2008 & Richardson 2011), when it rotate in a slower speed, the pre-stressed
status, dont reach yield strength. When increasing the speed to maximum allowance
rotational velocity, the pre-stressed status is reduced. To sum up, with the validation
in above section, the result of stresses of multi-rim flywheel in maximum rotations is
corrected.

CONCLUSION
After calculation of 3 cases, we have a table to sum it up:
One-rim
Case

Aluminum

One-rim T300

Alloy
Maximum velocity can
withstand (rad/s)
Maximum energy
capacity (J)
Mass = (Kg)

Two-rim
T300

1453.25

2066.096

5564.14

37954J

47091.2

341549.27

5.155524

3.294324

3.294324

Table 11. Results of 3 analyzed cases

With these results, the best material is composite, which is clearly having better
properties. The fact is the composite flywheel has higher energy capacity than
normal isotropic material, in this paper, 1.4 times. Furthermore, with the composite
flywheel, the weight is comparatively lighter, again in this paper, 1.7 times. The next
thing is we should consider another criterion of flywheel structure, which multi-rim is

25

better than single-rim. Based on the results, with the same materials, same volume, a
two-rim flywheel can store almost 8 times when compare with one-rim.
Although there was a limitation of this paper that the interference fit effect for
assembling 2 rims together wasnt conduct successfully, the model was validated
with multi approach for each of case. The report is able to conclude that multi-rim
composite flywheel have the best performance as predicted. So that is the optimized
choice for materials and structure of a flywheel using for energy storage systems.
In general, all of two objectives is achieved. Due to time and size limit of this
paper, the variation of each case was not examined, like changing some material,
dimension
. These will be next objectives in a future report. Thank you for reading.
REFERENCE
Chamis, C. C., & Kiraly, L. J. (1976). Rim-spoke composite flywheels: Stress and vibration
analysis.
Ha, S. et al. (2012). Design optimization and fabrication of a hybrid composite flywheel
rotor. Composite Structures,94(11), 3290-3299.
Flybrid Systems LLP (2010-09-10). "Flybrid Systems". Flybrid Systems.
Jost,

K.

(1994).

The

Patriot's

hybrid-electric

drivetrain. Automotive

Engineering;(United

States), 102(12).
Kamf, T. (2012). High speed flywheel design: Using advanced composite materials.
Abrahamsson, J., Hedlund, M., Kamf, T., & Bernhoff, H. (2014). High-speed kinetic energy buffer:
Optimization of composite shell and magnetic bearings.Industrial Electronics, IEEE Transactions

on, 61(6), 3012-3021.

Tokahashi, I., Itoh, Y., & Andoh, I. (1989, October). Development of a new uninterruptible power
supply using flywheel energy storage techniques. InIndustry Applications Society Annual
Meeting, 1989., Conference Record of the 1989 IEEE (pp. 711-716). IEEE.
Ricardo (2011) Breakthrough in Ricardo Kinergy second generation high-speed flywheel
technology, Accessed 22 May 2015 <http://www.ricardo.com />

Wen, S., & Jiang, S. (2012). Optimum design of hybrid composite multi-ring flywheel
rotor based on displacement method. Composites Science and Technology, 72(9),
982-988.

26

APPENDIX B: Orthotropic stresses calculation code MATLAB of Kamf (2012)

function [ Ten_c, Ten_r, intf ] = Ten_cr_calc(x,idat,n)
% Syntax:
% [Hoop stress, Radial stress, Radial displacement] =
% Ten_cr_calc(Input pressures, Geometrical data, Rotational speed)
% Calculates hoop and radial stresses of a spinning cylinder
% using theory chapter equations.
% Extract values from idat
m = idat(8,5);
ri = idat(1,1:1:m);
ro = idat(2,1:1:m);
Ec = idat(3,1:1:m);
Er = idat(4,1:1:m);
v_cr = idat(5,1:1:m);
p = idat(6,1:1:m);
y = idat(8,1);
G_rz = idat(9,1:1:m);
% Preallocating vectors
u = zeros(1,m);
C1 = zeros(1,m);
C2 = zeros(1,m);
Q = zeros(1,m);
Q2 = zeros(1,m);
r = linspace(min(ri), max(ro), y);
Ten_c = zeros(1,length(r));
Ten_r = zeros(1,length(r));
intf = zeros(1,length(r));
b = zeros(1,m);
v_rc = zeros(1,m);
v_rz = zeros(1,m);
v_zr = zeros(1,m);
v_zc = zeros(1,m);
v_cz = zeros(1,m);
% Useful relations
w = 2*pi*n/60; % rpm to rad/sec conversion
% Calculate Possions ratios for a assumed transversely isotropic material
for k = 1:1:m
v_rc(k) = v_cr(k)*Er(k)/Ec(k);
v_rz(k) = Er(k)/(2*G_rz(k))-1;
v_zr(k) = v_rz(k);
v_zc(k) = v_rc(k);
v_cz(k) = v_cr(k);
end
% Make Young's modulus quotas and b quota
for k = 1:1:m
u(k) = sqrt((Ec(k)/Er(k))*((1-v_rz(k)*v_zr(k))/(1-v_zc(k)*v_cz(k))));
b(k) = (v_cr(k)+v_cz(k)*v_zr(k))/(1-v_zc(k)*v_cz(k));
end
% Make C1 C2 and Q=K constants
for k = 1:1:m
Q(k) = (p(k)*w^2*(3+b(k))/(u(k)^2-9)); %last part of 2.20
Q2(k) = p(k)*w^2*(u(k)^2+3*b(k))/(u(k)^2-9); % last part of 2.22
C1(k) = (((ri(k)*ro(k))^(u(k)))/(ri(k)^(2*u(k))-ro(k)^(2*u(k))))...
*(Q(k)*(ri(k)^3*ro(k)^(u(k))-ro(k)^3*ri(k)^(u(k))));
C2(k) = (Q(k)*(ro(k)^(u(k)+3)-ri(k)^(u(k)+3)))/(ri(k)^(2*u(k))-ro(k)^(2*u(k)));
end
% Calculate stresses using equations from theory eqs
% and displacements using general Hooke's law
k = 1;
for i = 1:1:length(r)
k = min(k,m-1); %Limit k to prevent index overflow.
if r(i) > ri(k+1)
k = k+1;
end
k = max(k,1); %Limit k to prevent index overflow.
% Calculate stresses
Ten_r(i) = C1(k)*r(i)^(-1-u(k))+C2(k)*r(i)^(-1+u(k))+Q(k)*r(i)^2;
Ten_c(i) = u(k)*(C2(k)*r(i)^(-1+u(k))-C1(k)*r(i)^(-1-u(k)))+Q2(k)*r(i)^2;
%Calculate displacements
intf(i) = r(i)*(Ten_c(i)/Ec(k)-v_cr(k)*Ten_r(i)/Er(k));
end
end

27